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Abstract 

We  devise  a  novel  method  for  treating  bubbles  in  incompressible  flow  that  relies  on  the  conservative 
advection  of  bubble  mass  and  an  associated  equation  of  state  in  order  to  determine  pressure  boundary  con¬ 
ditions  inside  each  bubble.  We  show  that  executing  this  algorithm  in  a  traditional  manner  leads  to  stability 
issues  similar  to  those  seen  for  partitioned  methods  for  solid-fluid  coupling.  Therefore,  we  reformulate  the 
problem  monolithically.  This  is  accomplished  by  first  proposing  a  new  fully  monolithic  approach  to  coupling 
incompressible  flow  to  fully  nonlinear  compressible  flow  including  the  effects  of  shocks  and  rarefactions,  and 
then  subsequently  making  a  number  of  simplifying  assumptions  on  the  air  flow  removing  not  only  the  non- 
linearities  but  also  the  spatial  variations  of  both  the  density  and  the  pressure.  The  resulting  algorithm  is 
quite  robust,  has  been  shown  to  converge  to  known  solutions  for  test  problems,  and  has  been  shown  to  be 
quite  effective  on  more  realistic  problems  including  those  with  multiple  bubbles,  merging  and  pinching,  etc. 
Notably,  this  approach  departs  from  a  standard  two-phase  incompressible  flow  model  where  the  air  flow  pre¬ 
serves  its  volume  despite  potentially  large  forces  and  pressure  differentials  in  the  surrounding  incompressible 
fluid  that  should  change  its  volume.  Our  bubbles  readily  change  volume  according  to  an  isothermal  equation 
of  state. 


1.  Introduction 

Numerical  methods  for  simulating  two-phase  incompressible  flows  have  received  enormous  attention  pop¬ 
ularized  by  [3,  33,  32].  Quite  often,  the  air  phase  can  be  assumed  to  have  small  inertia  relative  to  the  water 
phase  and  need  not  be  simulated,  e.g.,  when  studying  the  impact  of  ocean  waves  on  a  ship,  which  gives  rise 
to  free  surface  flows.  However,  such  flows  ignore  effects  like  air  entrainment  and  are  unsuitable  for  modeling 
bubble  dynamics.  To  prevent  the  air  bubbles  from  collapsing  unnaturally,  methods  have  been  proposed 
for  computing  a  pressure  inside  these  bubbles  by  tracking  their  volume  over  time  and  using  an  equation 
of  state  [30].  But  volume  can  change  radically  if  a  bubble  rises  a  long  distance,  merges  with  other  nearby 
bubbles,  or  breaks  up  into  smaller  bubbles.  Since  the  total  mass  of  air  inside  water  is  conserved  over  time 
(ignoring  mass  exchange  across  the  interface  such  as  vaporization),  we  propose  to  track  the  bubble  mass 
instead  which  avoids  these  difficulties. 

We  devised  our  method  by  first  proposing  a  rather  straightforward  approach  based  on  mass  tracking 
in  Section  2.  This  approach  suffers  from  stability  issues  which  have  characteristics  similar  to  partitioned 
(as  opposed  to  monolithic)  methods  for  solid-fluid  coupling,  see  for  example  [13,  26,  12]  and  the  references 
therein.  These  issues  can  be  alleviated  using  outer  iterations  on  the  partitioned  solver,  as  discussed  in 
Section  2.4,  although  this  can  require  ten  or  more  Poisson  solves  per  time  step  and  is  thus  computationally 
expensive.  Instead,  we  take  a  more  traditional  monolithic  approach  for  the  air-water  problem  similar  to  the 
solid-fluid  coupling  in  [13,  26,  12]  as  motivated  by  [19].  We  begin  by  revisiting  the  partitioned  solver  for 
incompressible  and  compressible  flow  from  [5]  and  devise  a  monolithic  solver  using  the  ideas  from  [19]  to 
couple  together  incompressible  flow  with  fully  nonlinear  compressible  flow  including  shocks  and  rarefactions. 
The  results  of  this  method  are  shown  in  Section  3  for  both  the  gamma  gas  law  and  an  isothermal  equation 
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of  state.  Then  in  Section  4  we  simplify  this  approach  greatly  so  that  it  is  in  line  with  the  straightforward 
bubble  simulation  method  of  Section  2.  This  is  achieved  by  setting  both  the  bubble  density  and  the  bubble 
pressure  to  be  spatially  constant,  although  time- varying. 

In  summary,  our  method  for  simulating  bubbles  in  free  surface  incompressible  flows  is  as  follows.  The 
incompressible  flow  is  treated  in  a  standard  fashion  with  an  unconditionally  stable  advection  method  [8,  27], 
an  implicit  method  for  viscosity,  and  the  inclusion  of  surface  tension  using  jump  conditions  as  proposed  in  [22, 
18].  For  large  air  bodies  connected  to  the  atmosphere,  we  use  a  free  surface  condition  of  the  atmospheric 
pressure  as  a  Dirichlet  boundary  condition  as  is  done  in  standard  free  surface  incompressible  flow  solvers. 
For  bubbles,  we  define  mass  inside  each  bubble  and  conservatively  advect  that  mass  using  [21]  with  the 
air  velocity,  but  also  ensuring  that  the  mass  stays  confined  within  the  zero  level  set  using  ideas  proposed 
in  [20].  Density  variations  within  the  bubble  are  ignored  setting  the  bubble  density  equal  to  the  total 
bubble  mass  divided  by  the  total  bubble  volume  which  is  computed  using  a  flood  fill  algorithm.  The  bubble 
pressure  is  then  found  using  an  equation  of  state  p  =  Bp,  although  simply  applying  this  pressure  as  a 
Dirichlet  boundary  condition  on  the  incompressible  fluid  leads  to  instabilities  symptomatic  of  partitioned 
methods.  Instead,  we  use  the  pressure  evolution  equation  [10]  and  the  ideas  proposed  in  [19]  in  order  to 
formulate  a  fully  coupled  solve  between  a  single  degree  of  freedom  pressure  for  the  bubble  and  the  surrounding 
incompressible  flow.  This  provides  a  very  stable  monolithic  coupling  for  interactions  between  bubbles  and 
the  surrounding  incompressible  flow.  The  resulting  pressure  can  be  used  to  update  the  incompressible  flow  as 
usual,  whereas  the  air  velocities  are  updated  by  solving  a  second  Poisson  equation  with  Neumann  boundary 
conditions  forcing  the  bubble  velocities  to  match  the  motion  of  the  surrounding  incompressible  fluid  including 
both  expansion  and  contraction  allowing  the  bubble  to  change  volume  if  necessary.  Section  6  gives  a  more 
detailed  summary  of  our  method  referring  to  the  appropriate  equations  throughout  the  text  and  highlights 
its  efficacy  with  various  examples  -  including  the  ability  to  treat  multiple  bubbles  which  may  also  split  and 
merge. 


2.  A  partitioned  approach 

We  begin  with  a  straightforward  approach  that  conservatively  advects  the  bubble  mass  and  uses  the 
isothermal  equation  of  state  p  =  Bp  to  compute  a  pressure  inside  the  bubble  which  is  subsequently  used  as 
a  Dirichlet  boundary  condition  for  making  the  incompressible  velocities  divergence  free.  As  we  will  see,  this 
approach  lacks  stability  properties  because  the  current  pressure  in  the  bubble  is  not  a  good  predictor  of  what 
the  pressure  would  be  in  the  next  time  step  after  applying  the  incompressible  flow  pressure  to  generate  a 
new  velocity  field  that  changes  the  size  of  the  bubble.  If  one  couples  the  air  pressure  as  a  degree  of  freedom 
instead  of  as  a  Dirichlet  boundary  condition,  then  the  incompressible  Poisson  solver  can  better  react  to 
the  anticipated  changes  in  the  bubble  volume.  Thus,  after  presenting  our  straightforward  approach  in  this 
section  (which  we  call  partitioned),  we  present  a  monolithic  solver  in  Section  3. 


2.1.  Incompressible  flow 

The  incompressible  Navier- Stokes  equations  are  given  by 

_  _w  Vp  V  •  (flS/v) 

vt  +  (v-  V)v  H - =  - 

Pi  Pi 

V  •  v  =  0 


+  / 


(1) 

(2) 


where  pi  is  the  density,  v  is  the  velocity,  p  is  the  coefficient  of  viscosity  and  /  is  the  net  body  force  acting 
on  the  incompressible  fluid.  These  equations  are  discretized  on  a  MAC  grid  using  the  projection  method  [6], 
where  we  first  explicitly  update 
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and  then  solve  for  the  pressure  via 


Vp  V  •  v* 
Pi  A  t 

in  order  to  update  the  intermediate  velocity  v*  as  follows 

v n+1  —  v*  Vp 

At  +  pj 


(4) 

(5) 


We  use  the  level  set  method  [24]  to  track  the  interface  between  the  bubbles  and  the  incompressible  fluid. 
Before  updating  the  incompressible  velocities  through  advection,  they  are  extrapolated  across  the  interface 
in  order  to  define  ghost  node  values.  This  could  be  accomplished  using  constant  extrapolation  normal  to 
the  interface  by  solving  the  equation  IT  +  N  •  V/  =  0,  in  fictitious  time  r  for  each  component  /  of  v. 
We  instead  compute  the  steady  state  solution  using  the  fast  extension  method  of  [1].  The  incompressible 
velocities  are  then  advected  using  semi-Lagrangian  advection  which  can  be  made  second  order  accurate 
using  a  MacCormack-style  method  as  in  [27].  The  level  set  function  </>  is  advected  using  the  particle  level 
set  method  of  [7]  and  the  semi-Lagrangian  advection  scheme  of  [8].  To  keep  the  level  set  a  signed  distance 
function  we  use  the  modified  fast  marching  method  proposed  in  [23] . 

The  treatment  of  viscosity  for  multiphase  incompressible  flow  with  appropriate  jump  conditions  at  the 
interface  is  discussed  in  [18].  However,  viscosity  is  solved  for  explicitly  in  [18]  which  has  a  severe  time  step 
restriction  of  At  oc  0( Ax2).  In  order  to  take  large  time  steps,  we  consider  an  implicit  treatment  of  viscosity. 
As  discussed  in  [17],  if  all  the  jump  conditions  are  treated  implicitly  then  the  equations  for  all  components  of 
the  velocity  are  coupled  together.  Although  one  could  take  an  approach  similar  to  [25]  for  spatially  varying 
viscosity  where  the  coupling  terms  are  treated  explicitly  and  other  terms  are  treated  implicitly  in  order  to 
get  decoupling  of  various  components,  there  can  be  some  time  step  restrictions  based  on  the  jump  conditions. 
For  the  simulation  of  bubbles,  we  assume  that  the  dynamics  inside  the  bubbles  contain  little  momentum, 
hence,  they  cannot  absorb  any  viscous  momentum  from  the  liquid.  Thus,  we  enforce  Neumann  boundary 
conditions  at  the  interface  that  the  derivative  of  each  component  of  the  incompressible  velocity  is  zero.  Thus, 
the  jump  in  pressure  due  to  viscosity  is  also  zero  since  the  normal  component  of  the  viscous  stress  vanishes 
across  the  interface.  Finally,  as  we  assume  constant  viscosity  in  the  incompressible  fluid,  the  equations  for 
the  different  components  of  the  incompressible  velocity  decouple  as  well.  In  two  spatial  dimensions,  equation 
(3)  can  be  written  component-wise  as 
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The  advection  and  external  forces  can  be  applied  first  to  obtain  Vi  and  h2  followed  by  a  viscous  solve  of  the 
form 


v*  —  v\  V  •  (/iVt’i) 
At  P! 

V2  ~  V2  =  V  •  QV^) 

At  P! 


(8) 

(9) 


Since  Vi  and  h2  are  not  divergence  free  and  the  viscous  update  equations  have  been  derived  assuming 
the  divergence  free  condition,  v\  and  f)2  are  sometimes  first  projected  to  be  divergence  free  before  applying 
viscosity.  However,  note  that  the  pressure  projection  is  not  idempotent  in  the  presence  of  a  non-zero  pressure 
gradient  across  the  incompressible  fluid.  In  this  case  we  do  not  project  v\  and  h2  to  be  divergence  free  before 
the  viscous  update.  Note  that  the  advection  terms  are  computed  in  a  thin  band  of  ghost  cells  so  that  there 
are  adequate  values  when  the  interface  moves,  however,  the  viscous  terms  can  only  be  updated  interior  to 
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the  level  set  due  to  the  need  to  prescribe  interface  boundary  conditions.  Therefore,  the  level  set  must  be 
moved  to  its  new  time  tn+1  location  before  applying  the  viscous  update. 

We  use  the  second  order  cut-cell  method  of  [11]  for  computing  the  pressure  in  equation  (4)  to  make  the 
incompressible  velocities  divergence  free,  where  the  pressure  inside  the  bubbles  and  the  outside  air  are  used 
as  Dirichlet  boundary  conditions.  In  the  presence  of  surface  tension,  the  term  an  is  added  to  the  Dirichlet 
pressure,  where  a  is  the  coefficient  of  surface  tension  and  n  is  the  curvature  of  the  interface,  computed  using 
the  level  set  method  [24].  See  also  Section  5. 

The  incompressible  time  step  A tj  is  computed  by  enforcing  the  following  inequality  at  every  cell  center, 
as  described  in  [18], 


A  tj 


Ccfl  +  y/ Ccfl  +  4Scfl  +  4i^fl  j  <  1 


Here,  Cc fi  accounts  for  the  convection  terms  where  v\  and  have  been  averaged  to  the  cell  center, 


n  M  ,  H 

cfl  Ax  Ay 


San  accounts  for  the  surface  tension  forces, 

5cfi  = 


p/(min{Aa;,  Ay})2 


(10) 


(11) 


(12) 


and  Fcf|  accounts  for  the  body  forces  / 


Fcf,= 


'IM  +  IM 

Ax  Ay 


Note  that  the  time  step  restrictions  due  to  viscosity  are  not  present  as  it  is  treated  implicitly. 


(13) 


2.2.  The  oscillating  bubble  problem 

We  consider  a  model  oscillating  bubble  problem  in  one  spatial  dimension  as  shown  in  Figure  1,  where  an 
air  bubble  of  radius  r°  =  .1  m  with  initial  density  p°  =  1.1  kg/m3  is  inside  a  water  “sphere”  of  radius  =  .4 
m.  The  computational  domain  is  [0  m,  1  m]  which  gives  .1  m  of  free  air  on  each  side  of  the  water  region. 
Figure  2  shows  the  problem  in  two  spatial  dimensions  where  an  air  bubble  of  radius  r°  =  .1  m  with  initial 
density  p°  =  1.1  kg/m3  is  inside  a  water  sphere  of  radius  =  .4  m.  The  computational  domain  is  [0  m,  1 
m]  x  [0  m,  1  m].  The  setup  for  the  problem  in  three  spatial  dimensions  is  defined  similarly.  For  simplicity, 
there  is  no  gravity,  surface  tension  or  viscosity  in  the  system.  Since  the  bubble  is  slightly  compressed  with 
density  p°  =  1.1  kg/m3,  there  will  be  a  larger  pressure  p°  inside  the  bubble  than  in  the  ambient  air  which 
is  taken  to  be  a  free  surface  condition  of  pa tm  =  101,325  Pa  and  therefore,  the  bubble  will  start  to  expand, 
subsequently  vibrating  back  and  forth.  The  appendices  derive  a  second  order  ODE  given  by  equations  (63), 
(70)  and  (77).  In  all  three  equations,  we  take  the  standard  approach  of  solving  for  R(t),  rewriting  the  second 
order  equation  as  a  first  order  system,  subsequently  integrating  in  time  using  third  order  accurate  TVD 
Runge-Kutta,  and  refining  the  time  step  until  the  solutions  converge  to  obtain  data  that  we  use  for  the 
“exact”  solutions  when  these  equations  are  considered. 

For  this  problem,  we  modify  the  time  step  restriction  of  Section  2.1  to  account  for  velocities  near  zero 
when  the  bubble  volume  is  at  an  extrema,  i.e. ,  when  it  has  maximum  or  minimum  volume.  To  prevent  the 
time  step  from  becoming  excessively  large  in  these  cases,  we  add  a  term  to  At/  that  estimates  the  change  in 
velocity  over  a  time  step,  similar  in  spirit  to  what  was  done  in  [18]  for  body  forces  and  [19]  for  compressible 
flow  -  see  also  Section  3.4  for  a  quick  summary.  Essentially  what  is  needed  is  an  estimate  for  |Vp|  which  will 
influence  the  velocity.  We  do  this  by  computing  lx  and  ly  as  the  minimum  thicknesses  of  the  water  region 
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Air  Water  Bubble  Water  Air 


Figure  1:  Setup  for  the  oscillating  bubble  problem  in  one  spatial  dimension. 


in  the  x  and  y  directions,  and  approximate  \px\  as  |patm  ~  pn \/lx  and  \py\  as  |patm  —  pn\/ly,  where  pn  is  the 
pressure  inside  the  bubble  at  time  tn.  Then  we  write 


Figure  2:  Setup  for  the  oscillating  bubble  problem  in  two  spatial  dimensions. 


2.3.  Treating  bubbles 

Initially,  each  bubble  is  assigned  a  density  (or  mass),  and  the  density  is  advected  using  the  unconditionally 
stable  conservative  semi-Lagrangian  scheme  of  [21].  This  scheme  is  especially  effective  in  keeping  track  of 
small  bubbles,  since  the  level  set  loses  volume  over  time  and  cannot  keep  track  of  sub-grid  level  details. 
Although  one  could  advect  the  bubble  density  using  velocities  extrapolated  from  the  incompressible  flow, 
we  use  a  more  accurate  approach  where  air  velocities  are  constructed  and  maintained  in  a  separate  velocity 
field  u.  The  velocity  used  for  level  set  advection  is  a  hybrid  between  the  incompressible  flow  velocity  and  the 
air  velocities.  Even  so,  numerical  smearing  and  other  errors  will  cause  the  location  of  the  zero  level  set  and 
the  location  of  the  non-zero  bubble  densities  to  drift  apart  over  time.  We  address  this  as  follows.  First  we 
compute  the  total  mass  that  belongs  to  a  bubble  as  the  sum  of  all  the  mass  inside  the  bubble  and  all  the  mass 
closest  to  that  bubble.  Then  we  use  a  flood  fill  algorithm  on  that  bubble  to  identify  all  grid  cells  belonging 
to  that  connected  component.  The  volume  of  this  connected  component  is  carefully  computed  using  a 
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piecewise  linear  reconstruction  of  the  level  set  as  outlined  in  [20].  The  mass  is  then  uniformly  redistributed 
inside  the  bubble  to  obtain  a  spatially  constant  density  pt>-  The  Dirichlet  boundary  condition  used  for  the 
incompressible  Poisson  solve  given  in  equation  (4)  is  computed  using  the  equation  of  state  p  =  Bp  which 
simplifies  to  p  =  BM/V ,  where  M  and  V  are  the  total  mass  and  total  volume  of  a  bubble  respectively.  Here, 
B  is  taken  to  be  pat m  m3/kg,  so  that  a  density  of  1  kg/m3  yields  a  pressure  of  pat m  =  101,  325  Pa. 

Air  velocities  are  treated  in  a  manner  similar  to  the  incompressible  flow  velocities.  Ghost  node  values 
are  defined  using  the  fast  extension  method  of  [1]  exactly  as  is  done  for  the  incompressible  velocities,  and 
then  the  air  velocities  are  advected  using  the  MacCormack  method  [27].  Since  the  air  is  not  strictly  volume 
preserving  as  bubbles  can  expand  and  contract,  we  do  not  make  the  air  flow  incompressible  as  would  be 
done  for  a  truly  two-phase  incompressible  flow.  Instead,  we  solve  a  modified  form  of  equations  (4)  and  (5) 
for  each  bubble 


v  •  Vp  =  V  •  fT  -  V  •  ^+1  (15) 

fln+1  —  v*  +  Vp  =  0  (16) 

where  p  =  Atp/pb,  Pb  is  the  spatially  constant  air  density  inside  the  bubble  (which  could  be  different  per 
bubble),  and  u*  is  the  post-advected  air  velocity.  Equations  (4)  and  (5)  are  solved  first  for  the  full  water 
volume  using  Dirichlet  boundary  conditions,  after  which  the  resulting  water  velocities  surrounding  the  bubble 
are  used  as  Neumann  boundary  conditions  to  solve  equation  (15).  Since  the  integral  of  the  incompressible 
velocity  around  the  surface  of  a  bubble  may  not  be  zero,  we  compute  the  net  divergence  for  the  boundary 
of  the  bubble  summing  all  these  velocities  divided  by  the  number  of  cells  in  the  bubble,  and  use  that  value 
for  V  •  f2n+1.  This  allows  bubbles  to  expand  and  contract  and  otherwise  change  volume  as  they  follow  and 
are  enslaved  by  the  surrounding  incompressible  flow  which  has  much  higher  momentum.  Also  note  that  the 
Poisson  matrix  resulting  from  equation  (15)  has  a  rank- deficiency  of  1  due  to  the  full  Neumann  boundary 
conditions  and  although  the  addition  of  V  •  fln+1  guarantees  that  the  right  hand  side  is  in  the  range  of  the 
Poisson  matrix,  one  still  needs  to  take  care  to  compute  the  minimum  norm  solution  during  the  conjugate 
gradient  solve. 

In  Figures  3  and  4  we  present  numerical  profiles  generated  by  the  partitioned  scheme.  Consider  the  setup 
shown  in  Figure  1  and  assume  the  bubble  is  initially  at  rest,  i.e.,  the  initial  air  density  inside  it  is  p°  =  1 
kg/m3.  Analytically,  the  bubble  should  just  stay  at  rest  and  the  bubble  volume  should  remain  constant  over 
time.  We  refer  to  this  problem  as  the  stationary  bubble  problem.  Figure  3  shows  the  bubble  volume  profiles 
over  time  with  the  proposed  partitioned  scheme.  Note  that  although  the  initial  solution  computed  by  this 
scheme  is  indeed  constant,  it  quickly  goes  unstable  at  lower  grid  resolutions  because  of  the  forward  Euler 
characteristic  of  the  partitioned  scheme.  For  higher  grid  resolutions,  the  nonphysical  growth  is  slower  as 
expected.  For  the  one  dimensional  oscillating  bubble  problem,  Figure  4  shows  the  resulting  volume  profiles 
over  time  under  grid  refinement.  Again,  note  that  the  partitioned  scheme  shows  instability  at  lower  grid 
resolutions  while  converging  under  grid  refinement.  The  results  for  the  two  dimensional  simulations  are 
similar. 

2.4-  An  iterative  approach 

Although  the  partitioned  scheme  proposed  in  Section  2  converges  to  correct  analytical  solutions  under 
grid  refinement,  it  shows  instability  at  lower  grid  resolutions  because  of  its  forward  Euler  characteristic.  The 
most  logical  next  step  is  to  try  using  TVD  Runge-Kutta  methods  for  greater  stability  [28].  We  approach 
this  by  applying  second  order  accurate  TVD  Runge-Kutta  on  the  bubble  pressure  pb  which  is  used  as  a 
Dirichlet  boundary  condition  in  equation  (4).  Specifically,  we  take  two  full  steps  of  our  method  to  compute 
the  spatially  constant  density  p^+2  inside  the  bubble  and  use  the  average  (p^  +  p^+2)/2  for  computing  the 
Dirichlet  pressure  p b  =  B(p^  +  p^+2)/ 2.  Subsequently,  we  rewind  the  simulation  to  the  beginning  of  the 
time  step  and  use  this  Dirichlet  pressure  boundary  condition  in  the  incompressible  flow  solve  in  equations  (4) 
and  (5)  in  order  to  obtain  the  divergence  free  incompressible  velocity  field.  Note  that  intermediate  substeps 
can  result  in  velocities  that  dictate  a  smaller  step  size  than  that  chosen  at  the  beginning  of  the  time  step. 
Ignoring  this  can  result  in  inaccurate  solutions,  so  if  this  occurs  we  revert  to  the  beginning  of  the  time  step 
and  start  over  using  smaller  step  size.  Figure  5  shows  a  comparison  between  the  forward  Euler  scheme, 
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Forward  Euler:  CFL  =  .5,  initial  density  =  1. 


Figure  3:  Numerical  profiles  generated  by  the  partitioned  scheme  in  one  spatial  dimension  for  bubble  volume  over  time  for 
the  stationary  bubble  problem.  Note  the  forward  Euler  characteristic  of  the  partitioned  scheme  which  causes  instability  in  the 
computed  solutions  at  lower  grid  resolutions.  For  higher  grid  resolutions,  the  nonphysical  growth  is  slower  as  expected. 


second  order  accurate  TVD  Runge-Kutta  on  the  Dirichlet  pressure  pb  without  modifying  the  time  step, 
and  the  modified  second  order  accurate  TVD  Runge-Kutta  which  reverts  the  simulation  to  the  beginning 
of  the  time  step  if  the  CFL  condition  is  violated.  Note  the  intermittent  spikes  generated  by  the  standard 
second  order  accurate  TVD  Runge-Kutta  scheme,  which  are  not  present  in  the  modified  rewinding  version 
as  it  obeys  the  CFL  time  step  restriction  even  at  intermediate  steps.  Also  note  that  the  modified  second 
order  accurate  TVD  Runge-Kutta  scheme  does  a  much  better  job  at  tracking  the  constant  solution  than  the 
forward  Euler  version  of  the  scheme. 

The  enhanced  stability  shown  by  the  use  of  a  second  order  accurate  TVD  Runge-Kutta  scheme  on  the 
Dirichlet  pressure  boundary  condition  pb  used  for  the  bubble  in  solving  equations  (4)  and  (5)  motivates 
consideration  of  a  technique  similar  to  [19]  where  our  goal  is  to  replace  the  p  =  BM/V  Dirichlet  boundary 
condition  with  the  pressure  the  bubble  would  have  in  the  next  time  step  after  solving  for  the  incompressible 
velocities  and  advecting  the  bubble  forward  in  time.  Consider  the  pressure  evolution  equation  [10], 

pt  +  u  •  Vp  =  — pc2  V  •  u  (17) 

Since  we  assume  that  the  air  density  inside  the  bubble  is  spatially  constant,  it  follows  from  p  =  Bp  that 
Vp  =  0  and  thus 


pt  =  -pc2V  ■  u 


(18) 
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Forward  Euler:  CFL  =  .125,  initial  density  =  1.1 
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Figure  4:  Numerical  profiles  generated  by  the  partitioned  scheme  in  one  spatial  dimension  for  bubble  volume  over  time  for  the 
oscillating  bubble  problem.  Note  that  the  partitioned  scheme  shows  instability  at  lower  at  lower  grid  resolutions  because  of  its 
forward  Euler  characteristic,  although  converges  to  the  “exact”  solution  under  grid  refinement. 


For  a  gas  governed  by  the  equation  of  state  p  =  Bp,  the  sound  speed  c  is  defined  as 


PP  +  f  =  VB 

implying  that  it  is  constant  in  time  and  space.  Discretizing  equation  (18)  in  time  gives 

pn+1  =  pn  -  A tpn+1BW  •  u 


(19) 


(20) 


where  we  set  p  to  time  £n+1.  As  mass  is  conserved  over  time  Mn+1  =  Mn  =  M.  Using  p  =  Bp ,  pn  =  M /Vn 
and  pn+1  =  M/Un+1,  we  obtain 


BM  BM 


yn+l  yn 


BM 

’  yn+l 


A  tV  •  u 


(21) 


or 


AU  =  UnA tV  •  u  (22) 

where  AU  =  Un+1  —  Vn .  Using  the  divergence  theorem,  UnV  •  u  inside  the  bubble  is  equivalent  to  §  u-ndS. 
Letting  u  be  the  average  normal  velocity  on  the  boundary  of  the  bubble  and  V  be  the  perimeter  of  the 
bubble  allows  us  to  write 


AU  =  A  tVu 


(23) 
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Various  Methods:  CFL  =  .5,  initial  density  =  1. 


Figure  5:  A  comparison  between  numerical  profiles  for  the  one  dimensional  stationary  bubble  problem  generated  using  the 
forward  Euler  scheme  (red),  second  order  TVD  Runge-Kutta  (green),  modified  second  order  TVD  Runge-Kutta  which  takes 
time  step  restrictions  imposed  by  intermediate  velocities  into  account  (blue),  and  Brent’s  method  (magenta). 


Note  that  we  use  a  piecewise  linear  reconstruction  of  the  level  set  as  described  in  [20]  for  computing  V  and 

u. 

The  initial  guess  for  the  iterative  solver  sets  p  =  BM/V  as  a  Dirichlet  boundary  condition  for  projecting 
the  incompressible  velocities.  These  velocities  are  then  used  to  guess  the  bubble’s  volume  Vn+1  using 
equation  (23).  Since  mass  is  constant,  this  predicts  a  new  pn+1  =  Bpn+1  =  BM  /Vn+1 .  This  pressure  can 
again  be  set  as  a  Dirichlet  boundary  condition  to  project  the  incompressible  velocities  and  improve  the  guess 
for  pn+1.  Through  the  iterative  solver,  we  are  looking  for  a  fixed  point  for  this  “function”,  i.e.,  p  =  f(p)  or 
a  root  of  g(p)  =  f(p)  —p.  Note  that  if  the  input  pressure  is  too  large,  the  bubble  expands  and  the  predicted 
pressure  drops,  and  similarly  if  the  input  pressure  is  too  small,  the  bubble  contracts  and  the  predicted 
pressure  increases.  This  allows  us  to  place  bounds  on  the  solution.  Basically,  we  start  with  pn  =  BM/Vn 
and  if  f(p)  is  bigger,  the  initial  guess  for  p  is  the  left  end-point  of  our  interval.  Otherwise,  if  f(p)  is  smaller, 
this  is  taken  as  the  right  end-point  of  the  interval.  As  the  iteration  proceeds,  we  eventually  identify  both 
left  and  right  end-points.  Once  we  find  a  bounding  interval  we  use  Brent’s  method  [2]  to  find  the  root. 

Figures  6(a)  and  6(b)  show  the  numerical  profiles  generated  by  the  iterative  solver  for  the  bubble  volume 
over  time  with  increasing  grid  resolutions  in  one  and  two  spatial  dimensions  for  the  stationary  bubble 
problem.  Note  the  improved  stability  achieved  by  the  iterative  solver.  For  the  sake  of  comparison,  the 
iterative  solver  is  also  shown  labeled  as  Brent’s  method  in  Figure  5.  Next  consider  the  oscillating  bubble 
problem.  Figures  6(c)  and  6(d)  show  the  corresponding  profiles  for  the  bubble  volume  over  time  under  grid 
refinement.  Comparing  Figures  4  and  6,  note  that  the  explicit  method  diverges  for  low  grid  resolutions 
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Brent's  Method:  CFL  =  .125,  initial  density  =  1. 


Brent's  Method:  CFL  =  .5,  initial  density  =  1. 
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Brent's  Method:  CFL=.125,  initial  density=l.l 
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Figure  6:  Numerical  profiles  generated  by  the  iterative  scheme  for  bubble  volume  over  time  under  grid  refinement  for,  (a)  the 
stationary  bubble  problem  in  ID,  (b)  the  stationary  bubble  problem  in  2D,  (c)  the  oscillating  bubble  problem  in  ID,  (d)  the 
oscillating  bubble  problem  in  2D.  Note  the  improved  stability  achieved  by  the  iterative  method  as  compared  to  the  forward 
Euler  scheme. 


and  eventually  converges  to  the  “exact”  solutions  as  the  grid  is  refined.  In  contrast,  the  implicit  solution 
overly  damps  the  phenomena  on  coarse  grids  while  still  converging  to  the  “exact”  solutions  as  the  grid  is 
refined.  Obviously  one  would  prefer  the  implicit  approach  over  the  explicit  one  so  that  although  unresolved 
phenomena  are  not  accurately  resolved,  i.e.,  the  solution  is  overdamped,  they  at  least  do  not  explode  and 
corrupt  the  entire  solution. 


3.  Coupling  compressible  and  incompressible  flow 

Motivated  by  the  method  in  [19]  which  proposes  an  elliptic  solver  for  pressure  for  compressible  flow, 
our  goal  is  to  develop  a  fully  monolithic  solver  which  solves  for  air  and  water  pressures  together.  However, 
before  doing  this  we  first  develop  a  full  monolithic  solver  that  couples  incompressible  flow  with  fully  non¬ 
linear  compressible  flow  including  contact  discontinuities  like  shocks  and  rarefactions  building  upon  the 
partitioned  solver  of  [5].  They  used  a  fully  explicit  scheme  for  the  compressible  fluid,  and  the  pressure  in 
the  compressible  region  was  used  as  a  Dirichlet  boundary  condition  when  solving  the  incompressible  Poisson 
equation  for  updating  the  incompressible  velocities.  Near  the  interface,  the  ghost  fluid  method  (GFM)  [9]  was 
used  to  treat  the  boundary  conditions  in  a  manner  that  admitted  sharp  discontinuities  while  still  allowing 
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for  smooth  discretizations  across  the  interface.  To  achieve  this,  the  interface  values  of  pressure  and  normal 
velocity  were  carefully  determined  noting  that,  although  these  variables  are  continuous,  they  may  possess 
kinks  across  the  interface.  Small  errors  in  the  normal  velocity  of  the  incompressible  fluid  create  small  errors 
in  its  divergence,  which  in  turn  can  lead  to  large  spurious  pressure  oscillations  in  the  incompressible  region. 
While  small  errors  in  the  velocity  of  the  compressible  fluid  cause  the  same  small  errors  in  its  density,  these 
have  little  effect  on  the  gas  since  the  gamma  gas  law  equation  of  state  is  rather  robust.  Again,  since  the 
incompressible  flow’s  pressure  response  is  rather  stiff,  one  can  expect  large  variations  in  the  incompressible 
pressure  near  the  interface  which  in  turn  can  lead  to  poor  predictions  of  the  interface  pressure.  While  these 
errors  in  the  interface  pressure  have  a  relatively  small  effect  on  the  heavier  incompressible  fluid,  they  can 
have  a  rather  large  effect  on  the  lighter  compressible  gas.  Conversely,  since  the  gamma  gas  law  equation  of 
state  is  rather  robust,  the  compressible  pressure  tends  to  be  smooth  near  the  interface  and  is  therefore  a 
good  candidate  for  the  interface  pressure.  In  view  of  these  statements,  [5]  proposed  using  the  incompressible 
region  to  determine  the  interface  normal  velocity  and  the  compressible  region  to  determine  the  interface 
pressure.  In  the  presence  of  surface  tension,  the  compressible  pressure  is  not  used  directly,  but  is  first 
modified  according  to  the  appropriate  \p\  =  ok  jump  condition.  Although  the  method  of  [5]  works  well,  it 
suffers  from  a  strict  time  step  restriction  because  the  sound  speed  in  the  compressible  fluid  dictates  the  size 
of  the  time  step.  In  addition,  the  method  of  [5]  utilizes  a  partitioned  coupling  approach  which  can  suffer 
from  stability  issues. 


3.1.  A  semi-implicit  formulation  for  compressible  flow 

Let  p  be  the  density  of  the  compressible  fluid,  u  its  velocity,  E  the  total  energy,  p  the  pressure,  and 
U  =  (p,  pf?,  E )  the  compressible  state  vector.  The  inviscid  compressible  Euler  equations  in  multiple  spatial 
dimensions  are  as  follows 

(  p  \  /  pu  \  (  0  \ 

Ut  +  V  •  F(U)  =  pu  +  V  •  pu®u+p  =  0  (24) 

\  E  jt  V  (E+P)*  )  W 

A  semi-implicit  formulation  for  solving  these  equations  was  recently  proposed  in  [19],  where  the  flux  vector 
F(U)  was  split  into  an  advection  part  and  a  nonadvection  part 


Fi(U) 


pu  \  f  0  \ 

pu®u  ,F2(U)  =  p 
Eu  )  \  pu  ) 


(25) 


The  advection  part  Fi(U)  was  integrated  explicitly  to  get  intermediate  values  p*,  (pu)*  and  E*,  and  since 
pressure  does  not  affect  the  continuity  equation,  it  follows  that  pn+1  =  p* .  The  original  method  of  [19]  used 
a  modified  ENO  scheme  for  the  explicit  update  to  avoid  Gibbs  phenomena.  Recently,  [12]  proposed  using 
the  standard  ENO  scheme  [29]  and  a  different  method  for  computing  the  post-ad vected  pressure  (pa,  see 
below)  which  avoids  this  issue.  We  follow  this  improved  approach.  The  nonadvection  momentum  and  energy 
updates  are 

(26) 
(27) 

Motivated  by  the  standard  incompressible  flow  formulation,  equation  (26)  is  divided  by  pn+1  to  obtain 


( pu)n+1  —  (pu)* 
At 

En+i  _  /.;• 

A t 


=  — Vp 


=  —V  •  (pu) 


and  its  divergence  is  taken  to  obtain 


un+1 


=  u*  -At 


Vpn+ 1 

pn+1 


v  •  Mn+1  =  V  •  u*  -  AtV  • 


/Vpra+1\ 

V  Pn+1  ) 


(28) 


(29) 
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Then  the  pressure  evolution  equation  (17)  is  semi-discretized  by  fixing  V  •  u  to  time  £n+1  through  the  time 
step  and  by  treating  the  advection  terms  explicitly.  Let  e  =  E  /  p  —  u-u/2  denote  the  internal  energy  per  unit 
mass,  then  the  advected  pressure  is  computed  as  pa  =  P*  =  p{p *,  e*)  using  the  equation  of  state.  Substituting 
p*  into  the  semi-discretized  form  of  equation  (17),  we  obtain 

pU+l  =P*_  At/?C2V  .  $n+i  (30) 

Eliminating  V  •  iP+1  by  combining  equations  (30)  and  (29)  and  rearranging  the  terms  gives 

pn+i  _  A t2pn(c2)nX7  ■  (^^rj  =  P*~  a tpn{c2)nv  ■  u*  (31) 

where  the  term  pc 2  has  been  fixed  to  the  time  tn  value.  By  composing  the  pn(c2)n  terms  into  a  diagonal 
matrix  P  =  [At2pn(c2)n]  and  discretizing  the  gradient  and  divergence  operators,  we  obtain  the  following 
system  of  equations 


[P-1  +  GT(pn+1)-1G\pn+1  =  P- V  +  GTuk  (32) 

where  G  is  the  discretized  gradient  operator,  —  GT  the  corresponding  discretized  divergence  operator.  The 
pressure  is  scaled  by  At,  i.e.,  p  =  pAt,  p  is  the  density  interpolated  to  cell  faces  and  u  denotes  a  density- 
weighted  averaged  face  velocity,  i.e., 


pn+l  _  P_ 

P+1/2  “ 


n+1 


Pit I 


(pu)*  +  (, pu)t 


i+ 1/2 


i+ 1 


n+1 


p?t% 


(33) 


Note  that  the  resulting  matrix  in  equation  (32)  has  an  identity  term  in  it,  which  allows  fast  solvers  like 
preconditioned  conjugate  gradient  (PCG)  to  converge  in  relatively  few  iterations.  After  solving  equation 
(32)  to  obtain  cell-centered  pressure  values,  they  are  applied  in  a  conservative  flux-based  manner  to  update 
the  intermediate  momentum  and  energy.  Face  pressures  are  computed  using  density- weighted  averaging,  i.e., 


hn+1 

P+1/2 


p?+1+p^i 


(34) 


and  face  velocities  are  computed  by  rewriting  equation  (28)  using  face-averaged  quantities  as  defined  above. 


K+l/2  =  *1+1/2  -  (+++11/2)-1++l/2P 


n+1  1—1/ 


-n+1 


(35) 


Here,  1/2  denotes  the  row  of  G  corresponding  to  face  i  +  1/2.  The  flux-based  implicit  update  then  takes 
the  form 


(pu)^1  =  0 on )? 


hn+1 
Pi+ 1/2 


-  hn+1 
Pi- 1/2 


Ax 


E: 


n+1 


=  Et 


-n+1  -n+1 

P+1/2  z+1/2 


-n+1  -n+1 

P-1/2  z-1/2 


Ax 


(36) 


3.2.  Explicit  coupling  step 

We  use  the  level  set  method  to  track  the  interface  between  the  compressible  and  incompressible  fluids. 
For  the  explicit  part  of  the  method,  the  interface  boundary  conditions  are  treated  as  described  in  [5].  Various 
quantities  need  to  be  extrapolated  across  the  interface  in  either  direction  to  define  ghost  node  values.  This  is 
accomplished  using  constant  extrapolation  normal  to  the  interface  by  solving  the  equation  IT  +  N  •  V/  =  0, 
in  fictitious  time  r  for  the  different  quantities  I  -  we  use  the  fast  extension  method  of  [1].  To  compute  the 
compressible  ghost  node  state  we  decompose  the  extrapolated  state  vector  into  entropy,  pressure  and  velocity, 
compute  the  cell  centered  incompressible  velocity,  and  replace  the  normal  component  of  the  compressible 
velocity  field  with  the  normal  component  of  the  incompressible  cell  centered  velocity  field  before  reassembling 
entropy,  pressure  and  velocity  to  obtain  a  ghost  state  for  the  compressible  flow.  The  explicit  update  for  the 
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compressible  fluid  consists  of  applying  only  the  advection  fluxes  from  equation  (25)  using  these  ghost  node 
values.  In  order  to  properly  handle  uncovered  cells  a  band  of  ghost  node  values  is  also  updated  in  time. 

The  incompressible  flow  update  proceeds  similarly  as  to  what  was  described  in  Section  2.1.  The  full 
incompressible  velocity  field  is  extrapolated  across  the  interface  into  the  compressible  region  to  obtain  ghost 
node  values,  the  advection  terms  are  updated  to  obtain  intermediate  velocities,  and  the  viscous  terms  are 
updated  noting  that  the  compressible  fluid  is  inviscid  and  thus  cannot  absorb  any  viscous  stress. 

Note  that  when  applying  viscosity  to  the  incompressible  fluid,  the  time  t  velocities  Vi  and  v 2  are  sometimes 
first  projected  to  be  divergence  free,  as  noted  in  Section  2.1  after  equations  (8)  and  (9).  In  our  monolithic 
approach,  the  divergence  free  projection  is  two-way  coupled  between  the  compressible  and  incompressible 
flow  as  given  in  equations  (43)  and  (44)  below  in  Section  3.3.  Thus,  our  strategy  is  to  apply  this  two-way 
coupled  projection  first  and  then  use  the  resulting  incompressible  state  to  add  the  effects  of  the  viscous  terms 
to  the  incompressible  fluid.  Interestingly,  note  that  this  first  coupled  projection  solve  essentially  gives  the 
answer  one  would  obtain  if  the  incompressible  flow  was  inviscid.  Thus  for  the  viscous  case,  we  essentially 
obtain  the  inviscid  solution  first,  rewind  the  compressible  state  vector  U  to  its  pre-projected  state,  apply 
the  viscous  update  to  the  incompressible  region,  and  then  apply  the  coupled  projection  once  again  to  obtain 
a  final  solution  which  includes  the  effects  of  viscosity  on  the  incompressible  fluid. 

3.3.  Implicit  coupling  step 

For  the  sake  of  exposition,  we  describe  the  implicit  step  in  one  spatial  dimension.  Multiple  spatial 
dimensions  are  handled  in  a  straightforward  manner  using  a  dimension-by-dimension  approach.  Consider 
the  situation  depicted  in  Figure  7.  Let  p\nt  denote  the  interface  pressure  and  9  =  \cj)(x2)\/ (\4>(x2)\  +  \(/>(xs)\) 
be  the  cell  fraction  between  the  interface  and  the  center  of  cell  2.  Discretizing  the  incompressible  flow  Poisson 
equation  (4)  for  cell  3,  we  obtain 


Figure  7:  A  global  Poisson  solve  for  coupling  together  compressible  and  incompressible  fluids.  Compressible  cells  are  shown 
with  red  borders,  incompressible  cells  are  shown  with  blue  borders.  The  interface  is  shown  in  green  and  the  shared  face  is 
colored  black. 


A  t/p?+1-pg+1  ^ 

Ax  y  pi  Ax  pi(  1  —  0)  Ax ) 


(37) 


Let  I  and  C  subscripts  represent  values  from  the  incompressible  and  compressible  sides  of  the  interface.  For 
inviscid  flow  Vp/p  is  continuous  across  the  interface  [18]  satisfying 


Vp/  _  Vpc 

Pi  Pc 


(38) 


Although  [18]  showed  this  for  two-phase  incompressible  flow,  one  can  still  derive  equation  (1)  for  compressible 
flow  from  the  equations  for  conservation  of  mass  and  momentum  and  so  it  also  holds  for  inviscid  compressible 
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flow  and  mixed  compressible-incompressible  flow  (as  long  as  the  strong  form  of  the  equations  hold  and 
derivatives  exist).  For  viscous  flows  things  can  be  more  complex,  however,  in  our  case  we  assume  the 
compressible  flow  is  inviscid  and  that  there  is  no  viscous  momentum  transfer  from  the  incompressible  to  the 
compressible  flow  as  mentioned  in  Section  2.1  (see  also  Section  3.2)  -  allowing  equation  (38)  to  still  hold. 
Approximating  equation  (38)  with  one-sided  differences,  we  obtain 

r).  _  ^+1  ^+1  _  r). 

Pint  V2  _  V_ 3 _ Pint  /qq\ 

pcOAx  pi(l  —  0)  Ax 


where  we  take  pc  =  P2+1  (one  might  also  conceivably  use  p^1).  Solving  equation  (39)  for  the  interface 
pressure  p-mt  gives 


Pint  — 


Opcp +1  +  (i  -  %jp?+1 

Opc  +  (1  -  0)pi 


(40) 


Writing 


p  =  0pc  +  (  1  -  0)pi 


(41) 


allows  us  to  write 


Pint  ~P2+1 


P3  +  1 


Pn2  +  1 


P3+1  “Pint 


pcOAx  pAx  (1  —  0)piAx 


(42) 


and  thus,  equation  (37)  becomes 


At  /p"+1-pg+1  p^+1-P2+1 

Ax  \  pi  Ax  pAx 


-V  -v* 


(43) 


Next  consider  the  compressible  Poisson  equation  in  cell  2.  Note  that  cell  3  has  a  valid  compressible  state 
vector  after  advection  because  we  also  update  a  band  of  ghost  cells  near  the  interface  (see  Section  3.2).  Let 
P2  jj1,  Put1  denote  the  interpolated  face  density  obtained  by  averaging  cell-centered  values.  Discretizing  the 
compressible  Poisson  equation  (32)  for  cell  2  and  using  equation  (42)  gives 


^+1  _  *  (Pl t  -  ^+1  _  ^+1  -  PTl  \=  Pl  _  y .  (44) 

p%(cl)nAt  Ax  \  pAx  p7^1  Ax  )  p2{cl)nAt 

Note  that  equations  (43)  and  (44)  together  form  a  symmetric  positive  definite  (SPD)  system,  allowing  for 
the  use  of  fast  Poisson  solvers  such  as  PCG. 

Besides  using  p  in  equations  (43)  and  (44)  to  enforce  the  balance  condition  of  equation  (38)  we  also 
desire  a  unique  velocity  at  the  shared  face  location  £2.5  (shown  in  black  in  Figure  7)  for  computing  the 
term  V  •  v*  in  equation  (43)  and  the  term  V  •  u*  in  equation  (44).  Since  we  have  separate  velocity  fields 
for  compressible  and  incompressible  fluids,  the  values  U2.5  an(^  ^2.5  shared  face  need  not  be  equal, 

even  though  theoretically  the  normal  component  of  the  velocity  field  is  supposed  to  be  continuous  across  the 
interface.  Thus,  at  the  shared  face  we  apply  a  force  A  to  the  incompressible  fluid  and  an  equal  and  opposite 
force  —A  to  the  compressible  fluid  so  that  the  resulting  time  t**  velocities  at  the  shared  face  are  equal,  i.e., 


~k~k 


•71+ 1  _ 


P2.5  ^2.5  =  P2.5  ^2.5 


AAt,  PivZX  =  piv 9  c  +  AA t 


(45) 


Setting  ^2*5  —  V2%:  and  solving  for  A  gives 


x2.5 


V2.5  — 


PIV  2.5 


'  P2%lu‘ 


2.5 


Pi  +  P2.51 


(46) 


as  our  velocity  at  the  shared  face  (note  that  since  the  incompressible  density  tends  to  be  much  bigger  than 
the  compressible  density,  one  could  also  use  the  incompressible  velocity  and  obtain  similar  results  -  we  have 
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tested  this  numerically).  Also  note  that  the  compressible  face  velocity  is  not  an  actual  degree  of  freedom  since 
the  degrees  of  freedom  for  compressible  flow  he  at  cell  centers.  Hence,  we  add  the  momentum  — AA t  to  both 
the  compressible  cell  2  and  the  compressible  ghost  cell  3,  i.e.,  ( pu )£*  =  (pu)^  —  XAt  and  ( pu )£*  =  (pu)\  —  \Aty 
so  that  the  first  equation  of  equation  (45)  is  satisfied.  This  is  equivalent  in  spirit  to  adding  —  AA t  to  cell 
2  and  then  re-extrapolating  into  ghost  cell  3.  Finally,  after  changing  the  compressible  state  vector  at  cell 
2,  the  velocity  u\  5  is  recomputed  to  make  it  consistent  with  the  new  state  in  cell  2.  In  multiple  spatial 
dimensions,  if  a  compressible  cell  borders  multiple  incompressible  cells  then  it  would  be  updated  multiple 
times,  leading  to  the  compressible  velocity  at  the  shared  face  not  matching  the  incompressible  velocity.  We 
use  the  updated  incompressible  velocity  at  the  shared  face  in  this  case  for  greater  accuracy. 

After  solving  the  coupled  projection  solve  of  equations  (43)  and  (44),  the  incompressible  velocity  at  the 
shared  face  is  updated  to  time  £n+1  via 


n+ 1  _ 


v2.5 


—  ^2.5 


P3  +  1 


-Pn2+1 


pAx 


(47) 


For  updating  the  compressible  momentum  at  cell  2,  we  compute  the  face  pressure  pJs  using  the  linear 
interpolant  through  p^1  and  pjnt,  i.e.,  p^1  =  p 2+1  +  (pint  —  P2+1)/^-  If  #  is  too  small  for  the  denominator, 
we  instead  compute  the  slope  using  equation  (42)  and  obtain  p^1  =  P2+1  +  P2+1(^3+1  —  P2+1)/(2p)-  The 
compressible  momentum  at  cell  2  is  then  updated  as 


(pu)n2+1  =  (pu)? 


Ax 


(48) 


noting  that  ( pu )£*  includes  the  momentum  update  AA t.  For  the  energy  update,  equation  (36)  is  still  used 
with  the  face  pressure  p^1  as  computed  above  and  the  face  velocity  u7^1  at  the  shared  face  computed  using 
&2*5  in  equation  (35).  Note  that  in  the  presence  of  surface  tension,  we  modify  the  Poisson  equations  for  both 
the  compressible  cell  2  and  the  incompressible  cell  3  by  taking  the  jump  ok  into  account,  as  described  in 
Section  2.1  (see  also  [18,  22]). 


3.4-  Time  step  restriction 

The  size  of  the  overall  time  step  is  computed  as  the  minimum  of  the  incompressible  and  the  compressible 
time  steps,  i.e., 


At  =  a  min{At/,  Ate} 


(49) 


where  a  denotes  the  CFL  number.  The  incompressible  time  step  At/  is  computed  as  described  in  Section  2.1. 
The  compressible  time  step  Ate  is  computed  as  described  in  [19].  In  order  to  prevent  Ate  from  becoming 
infinite  for  near  zero  velocities  un ,  the  term  Vp/p  is  added  which  estimates  the  change  in  velocity  at  the  end 
of  a  time  step.  Hence,  the  CFL  condition  becomes 


Atc 


+  J^A  tc 


PC 


Ax 


<  1 


(50) 


Equation  (50)  is  a  quadratic  in  Ate  with  two  solutions 


^/Mmax  +4^Aa;  <  <  -|un|max  +  ^/|un|^ax  +  4^Aa; 

2 \px\/pc  -  At°  -  2\Px\/pc 


(51) 


Note  that  the  lower  limit  in  equation  (51)  is  non-positive,  and  as  Ate  >  0,  only  the  upper  bound  needs  to 
be  enforced.  As  px  — >  0,  both  the  numerator  and  the  denominator  vanish  obtaining  the  typical  bound  of 
Ax/ 1 id1 1  max  which  is  problematic  when  |^n|max  is  small.  We  obtain  a  more  convenient  time  step  restriction 
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which  is  not  plagued  by  either  small  \px\  or  small  |^n|max  by  replacing  the  second  At  in  equation  (50)  with 
the  right  hand  bound  from  equation  (51)  to  obtain 


Atc 

2 


|  Un  |  max 

Ax 


+ 


(52) 


In  two  spatial  dimensions,  the  following  CFL  restriction  is  obtained 

^ C  (  \ul  |  max  1^2  Imax  If  l^llmax  j^£|max\  ^  \Px\  ^  \Py\ 

2  1  Ax  Ay  y  \  Ax  Ay  )  pc  Ax  pc^U 

3.5.  Numerical  Results 

We  used  the  gamma  gas  law  equation  of  state  p  =  (7  —  l)pe  and  an  outer  loop  of  third  order  TVD 
Runge-Kutta  in  all  these  examples,  noting  that  although  the  implicit  treatment  of  terms  inside  the  RK  loop 
generally  leads  to  a  loss  of  third  order  time  accuracy,  greater  stability  is  achieved.  All  examples  use  a  CFL 
number  of  .5. 


3.5.1.  One  dimensional  examples 

Consider  a  computational  domain  of  [0  m,  1  m].  The  domain  is  filled  with  a  compressible  gas  with 
p  =  1.226  kg/m3,  u  =  0  m/s,  and  p  =  105  Pa.  An  incompressible  droplet  of  length  .2  m  is  located  at 
the  center  of  the  domain  with  p  =  1000  kg/m3,  u  =  100  m/s,  and  p  =  105  Pa.  Since  the  incompressible 
droplet  is  moving  rightwards  in  an  initially  stationary  gas,  a  shock  wave  forms  in  the  gas  ahead  of  it  and 
a  rarefaction  wave  forms  in  the  gas  behind  it.  Figure  8  shows  the  density,  velocity,  pressure  and  time  step 
profiles,  along  with  the  fully  explicit  method  of  [5]  run  on  the  finest  grid  of  resolution  12800  for  the  sake  of 
comparison.  (Note  that  our  implementation  of  [5]  gave  a  different  result  for  this  example,  although  agrees 
with  [5]  for  all  other  examples  that  they  ran  in  both  one  and  two  spatial  dimensions,  leading  us  to  believe 
that  there  is  probably  a  typo  in  the  description  of  this  example  in  [5].)  Note  that  our  method  converges  to 
the  highly  refined  explicit  result  as  the  grid  is  refined.  Figure  9  shows  similar  results  when  the  density  of 
the  incompressible  droplet  is  p  =  10  kg/m3.  The  compressible  gas  slows  down  the  lighter  droplet  faster  and 
as  a  result  secondary  rarefaction  waves  stretch  between  the  droplet  and  the  shock  and  rarefaction  waves. 

Note  that  our  method  can  take  a  time  step  that  is  four  times  larger,  however,  the  cost  of  each  time 
step  is  slightly  larger  because  the  compressible  degrees  of  freedom  have  been  added  to  the  incompressible 
Poisson  solver.  The  overall  speedup  in  wall  clock  time  will  generally  depend  on  the  ratio  of  increased  cost 
per  time  step  as  compared  to  the  decrease  in  the  number  of  time  steps  -  in  this  particular  example,  the  code 
was  approximately  three  times  faster.  For  some  problems  the  speedups  can  be  significantly  larger  especially 
when  one  cares  about  phenomena  that  occur  after  the  compressible  flow  relaxes  to  smaller  velocities  -  in  this 
case  the  time  steps  could  increase  by  several  orders  of  magnitude. 

Consider  a  computational  domain  of  [0  m,  1  m]  filled  with  a  compressible  fluid  with  density  p  =  1.58317 
kg/m3,  velocity  u  =  0  m/s,  and  pressure  p  =  98066.5  Pa.  An  incompressible  droplet  of  length  .2  m  is  initially 
located  at  the  center  of  the  domain  with  p  =  1000  kg/m3,  u  =  0  m/s,  and  p  =  98066.5  Pa.  A  shock  wave  is 
initially  located  at  x  =  .1  m  with  a  post-shock  state  of  p  =  2.124  kg/m3,  u  =  89.981  m/s,  and  p  =  148407.3 
Pa  to  the  left  of  x  =  .1  m.  The  shock  wave  travels  to  the  right  impinging  on  the  incompressible  droplet, 
causing  both  reflected  and  transmitted  waves  as  shown  in  Figure  10  at  t  =  1.75  x  10-3  s.  Note  that  the 
transmitted  wave  is  too  weak  to  be  seen  in  this  example,  although  it  can  be  clearly  seen  in  Figure  11  where 
the  incompressible  droplet  has  density  p  =  10  kg/m3. 

3.5.2.  Two  dimensional  examples 

All  two-dimensional  examples  include  the  effects  of  viscosity  and  surface  tension  with  coefficients  p  = 
.001137  kg/ms  and  a  =  .0728  kg/s2.  These  effects  are  not  present  in  the  one  dimensional  examples  shown 
in  Section  3.5.1  because  the  incompressible  flow  has  constant  velocity  and  the  interface  has  no  curvature. 
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Moving  Heavy  Incompressible  Droplet:  Density 


Moving  Heavy  Incompressible  Droplet:  Velocity 


(a)  density 

Moving  Heavy  Incompressible  Droplet:  Pressure 


(b)  velocity 


Moving  Heavy  Incompressible  Droplet:  Time  Step 


time 


(c)  pressure 


(d)  time  step 


Figure  8:  Numerical  results  for  the  moving  incompressible  droplet  example,  where  a  droplet  of  density  1000  kg/m3  is  travelling 
to  the  right  in  an  initially  stationary  compressible  fluid  at  t  =  7.5  X  10-4  s. 
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Moving  Light  Incompressible  Droplet:  Density  Moving  Light  Incompressible  Droplet:  Velocity 


(a)  density 


(b)  velocity 


Moving  Light  Incompressible  Droplet:  Pressure  Moving  Light  Incompressible  Droplet:  Time  Step 


time 


(c)  pressure 


(d)  time  step 


Figure  9:  Numerical  results  for  the  moving  incompressible  droplet  example,  where  a  droplet  of  density  10  kg/m3  is  travelling 
to  the  right  in  an  initially  stationary  compressible  fluid  at  t  =  7.5  X  10-4  s. 
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Shock  Impinging  on  Heavy  Incompressible  Droplet:  Density 


Shock  Impinging  on  Heavy  Incompressible  Droplet:  Velocity 


(a)  density 


(b)  velocity 


Shock  Impinging  on  Heavy  Incompressible  Droplet:  Pressure  Shock  Impinging  on  Heavy  Incompressible  Droplet:  Time  Step 


time 


(c)  pressure 


(d)  time  step 


Figure  10:  Numerical  results  for  the  shock  impinging  on  a  heavy  incompressible  droplet  example  at  t  =  1.75  X  10-3  s,  where 
a  shock  wave  initially  located  at  x  =  .1  m  travels  to  the  right  impinging  on  an  incompressible  droplet  of  density  1000  kg/m3 
generating  both  reflected  and  transmitted  waves. 
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Shock  Impinging  on  Light  Incompressible  Droplet:  Density 


Shock  Impinging  on  Light  Incompressible  Droplet:  Velocity 


(a)  density 

Shock  Impinging  on  Light  Incompressible  Droplet:  Pressure 


(c)  pressure 


(b)  velocity 

Shock  Impinging  on  Light  Incompressible  Droplet:  Time  Step 


time 

(d)  time  step 


Figure  11:  Numerical  results  for  the  shock  impinging  on  a  light  incompressible  droplet  example  at  t  =  1.75  X  10-3  s,  where 
a  shock  wave  initially  located  at  x  =  .1  m  travels  to  the  right  impinging  on  an  incompressible  droplet  of  density  10  kg/m3 
generating  both  reflected  and  transmitted  waves. 
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Figure  12:  The  moving  incompressible  droplet  example  in  two  spatial  dimensions,  where  an  incompressible  droplet  of  density 
1000  kg/m3  is  travelling  to  the  right  in  an  initially  stationary  compressible  fluid,  (a)  50  equally  spaced  pressure  contours 
between  .75  X  105  Pa  and  1.5  X  105  Pa  on  a  1600  X  1600  grid  at  t  =  5  X  10-4  s,  (b)  pressure  contour  of  1.1  X  105  at  t  =  5  X  10-4 
s  under  grid  refinement  to  illustrate  convergence,  (c)  velocity  field  at  t  =  5  X  10-4  s  where  the  incompressible  velocities  are 
shown  in  blue,  and  compressible  velocities  are  shown  in  red,  and  (d)  the  zero  level  set  under  grid  refinement  at  t  =  2.5  X  10— 3 
s  as  compared  to  its  initial  location. 
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Figure  13:  The  moving  incompressible  droplet  example  in  two  spatial  dimensions,  where  an  incompressible  droplet  of  density 
10  kg/m3  is  travelling  to  the  right  in  an  initially  stationary  compressible  fluid,  (a)  pressure  contour  of  1.1  X  105  at  t  =  5  X  10-4 
s  under  grid  refinement  to  illustrate  convergence,  (b)  the  zero  level  set  under  grid  refinement  at  t  =  2.5  X  10-3  s  as  compared  to 
its  initial  location,  (c)  the  zero  level  set  on  a  grid  of  resolution  800  X  800  at  t  =  2.5  X  10-3  s  as  compared  to  its  initial  location 
using  the  fully  explicit  partitioned  solver  of  [5],  and  (d)  the  zero  level  set  on  a  grid  of  resolution  800  X  800  at  t  =  2.5  X  10-3  s 
as  compared  to  its  initial  location  using  our  monolithic  solver. 
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Consider  a  computational  domain  of  [0  m,  1  m]x[0  m,  1  m].  Similar  to  the  one  dimensional  case,  the 
domain  is  filled  with  a  compressible  gas  with  p  =  1.226  kg/m3,  u  =  v  =  0  m/s,  and  p  =  105  Pa.  An 
incompressible  droplet  of  radius  .2  m  is  located  at  the  center  of  the  domain  with  p  =  1000  kg/m3,  u  =  100 
m/s,  v  =  0  m/s,  and  p  =  105  Pa.  Since  the  compressible  gas  is  initially  stationary  and  the  droplet  is  moving 
rightwards,  a  shock  wave  forms  in  the  gas  in  front  of  it,  and  a  rarefaction  wave  forms  in  the  gas  behind 
it.  Figure  12(a)  shows  50  equally  spaced  pressure  contours  between  .75  x  105  Pa  and  1.5  x  105  Pa  on  a 
1600  x  1600  grid  at  t  =  5  x  10  4  s.  Figure  12(b)  shows  the  pressure  contour  of  1.1  x  105  at  various  grid 
resolutions  to  show  that  the  numerical  profiles  generated  using  our  method  converge  to  those  generated  using 
the  fully  explicit  method  of  [5]  under  grid  refinement.  The  velocity  field  is  shown  in  Figure  12(c)  where  the 
incompressible  velocities  are  shown  in  blue  and  the  compressible  ones  are  shown  in  red.  Figure  12(d)  shows 
the  initial  location  of  the  zero  level  set  as  compared  to  its  location  at  t  =  2.5  x  10-3  s.  Figure  13(a)-(d)  show 
the  results  for  the  case  when  the  incompressible  droplet  has  p  =  10  kg/m3.  Note  that  the  lighter  droplet 
undergoes  larger  deformation  and  also  slows  down  at  a  faster  rate.  Also  note  that,  as  observed  in  [5],  the 
computations  for  the  lighter  droplet  show  signs  of  Kelvin- Helmholtz  instability  as  is  apparent  by  the  wiggles 
in  the  interface  location  shown  in  Figure  13(b).  This  effect  is  less  apparent  on  coarser  grids  because  of  the 
artificial  damping  due  to  numerical  viscosity. 


Figure  14:  The  moving  incompressible  droplet  example  in  two  spatial  dimensions,  where  an  incompressible  droplet  of  density  10 
kg/m3  is  travelling  to  the  right  in  a  small  domain  in  an  initially  stationary  compressible  fluid,  (a)  one  dimensional  cross-section 
of  the  pressure  on  an  800  X  800  grid  at  t  =  5  X  10-9  s,  where  the  pressure  in  the  incompressible  region  is  shown  in  blue  and 
that  in  the  compressible  region  is  shown  in  red,  and  (b)  the  zero  level  set  under  grid  refinement  at  t  =  2.5  X  10— 8  s  as  compared 
to  its  initial  location. 

To  demonstrate  the  effects  of  surface  tension  and  viscosity,  we  also  shrunk  the  domain  to  [0  m,  lxlO-5 
m]x[0  m,  lxlO-5  m]  for  the  case  when  the  incompressible  droplet  has  density  p  =  10  kg/m3.  Figure 
14(a)  shows  a  one  dimensional  cross-section  of  the  pressure  at  t  =  5  x  10-9  s,  where  the  pressure  in  the 
incompressible  region  is  shown  in  blue  and  that  in  the  compressible  region  is  shown  in  red.  Note  the  jump 
in  pressure  across  the  interface  due  to  surface  tension  effects.  Figure  14(b)  shows  the  initial  location  of  the 
zero  level  set  as  compared  to  its  location  at  t  =  2.5  x  10-8  s.  Note  that  the  smaller  droplet  is  deformed  less 
and  has  a  more  spherical  shape  as  compared  to  the  larger  droplet  shown  in  Figure  13(d). 

Consider  a  computational  domain  of  [0  m,  1  m]x[0  m,  1  m].  Similar  to  the  one  dimensional  case,  the 
domain  is  filled  with  a  compressible  fluid  with  p  =  1.58317  kg/m3,  u  =  v  =  0  m/s,  and  p  =  98066.5  Pa.  An 
incompressible  droplet  with  initial  state  p  =  10  kg/m3,  u  =  v  =  0  m/s,  and  p  =  98066.5  Pa  with  radius  .2  m 
is  located  at  the  center  of  the  domain.  A  shock  wave  is  initially  located  at  x  =  .1  m  with  a  post-shock  state 
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Figure  15:  The  shock  impinging  on  an  incompressible  droplet  example  in  two  spatial  dimensions,  where  a  shock  wave  initially 
located  at  x  =  .1  m  travels  to  the  right  impinging  on  an  incompressible  droplet  with  p  =  10  kg/m3  generating  both  reflected  and 
transmitted  waves,  (a)  50  equally  spaced  contours  between  1  x  105  Pa  and  1.8  X  105  Pa  on  a  1600  X  1600  grid  at  t  =  1.25  X  10-3 
s,  (b)  pressure  contour  of  1.62  x  105  Pa  under  grid  refinement  at  t  =  1.25  X  10-3  s  to  illustrate  convergence,  (c)  velocity  field 
at  t  =  1.25  X  10— 3  s  where  the  incompressible  velocities  are  shown  in  blue,  and  the  compressible  velocities  are  shown  in  red, 
and  (d)  the  zero  level  set  under  grid  refinement  at  t  =  2.5  X  10-3  s  as  compared  to  its  initial  location. 
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of  p  =  2.124  kg/m3,  u  =  89.981  m/s,  v  =  0  m/s,  and  p  =  148407.3  Pa  to  the  left  of  x  =  .1  m.  The  shock 
wave  travels  to  the  right  impinging  on  the  incompressible  droplet,  generating  both  reflected  and  transmitted 
waves.  Figure  15(a)  shows  50  equally  spaced  pressure  contours  between  1  x  105  Pa  and  1.8  x  105  Pa  on 
a  1600  x  1600  grid  at  t  =  1.25  x  10-3  s.  Figure  15(b)  shows  the  pressure  contour  of  1.62  x  105  Pa  under 
grid  refinement  to  illustrate  convergence  to  those  generated  using  the  fully  explicit  method  of  [5].  Figure 
15(c)  shows  the  velocity  field,  where  the  incompressible  velocities  are  shown  in  blue,  and  the  compressible 
velocities  are  shown  in  red.  Figure  15(d)  shows  the  initial  location  of  the  zero  level  set  as  compared  to  its 
location  at  t  =  2.5  x  10-3  s  under  grid  refinment.  Note  that  the  computations  on  the  finer  grids  also  show 
signs  of  Kelvin- Helmholtz  instability. 


3. 6.  Constant  temperature  formulation 

We  make  an  isothermal  assumption,  where  the  equations  for  conservation  of  mass  and  momentum  form 
a  closed  system  and  the  equation  for  conservation  of  energy  decouples.  In  this  case,  the  gamma  gas  law 
equation  of  state  can  be  rewritten  as  p  =  Hp,  where  B  is  essentially  (7  —  l)e  which  we  set  equal  to  the 
atmospheric  pressure,  i.e.,  B  =  patm/Patm  m3/kg,  as  described  in  Section  2.  For  this  equation  of  state, 
substituting  equation  (19),  equation  (31)  becomes 

pn+ 1  -  A t2Bpn\7  •  (////)  =  Bp*  -  AtBpnX7  ■  u*  (54) 

where  p *  has  been  replaced  with  Bp*.  Note  that,  as  mentioned  in  Section  3.1,  p*  =  pn+1,  and  thus 
p*  =  Bp*  =  Bpn+l .  Furthermore,  if  pc 2  in  equation  (31)  is  taken  to  be  at  time  £n+1  instead  of  £n,  we  can 
replace  pn  with  pn+1  in  equation  (54)  and  divide  by  Bpn+1  to  obtain 

Sp7«#”+,-V  =  (55) 

which  is  the  continuous  analog  of  equation  (32)  for  this  equation  of  state. 


Figure  16:  Numerical  profiles  for  bubble  volume  over  time  for  the  oscillating  bubble  problem  in  one  spatial  dimension  with 
equation  of  state  p  =  Bp  generated  by  (a)  the  method  of  [5],  and  (b)  the  monolithic  solver. 


3.6.1.  Numerical  results 

Now  reconsider  the  examples  presented  in  Section  3.5.1  simulated  using  p  =  Bp  as  the  equation  of  state. 
We  use  the  same  ambient  conditions  for  density  and  velocity  in  all  the  examples,  noting  that  the  ambient 
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pressures  will  be  different  since  pressure  depends  on  density.  Also,  for  the  examples  with  prescribed  shocks 
we  choose  to  match  the  shock  speed  prescribing  a  post-shock  state  of  p  =  1.97705  kg/m3  and  u  =  70.4023 
m/s.  Figures  18-21  show  the  numerical  profiles  generated  using  equation  (55)  -  a  high  resolution  comparison 
to  the  fully  explicit  method  of  [5]  is  also  shown  in  the  results. 


Partitioned  Solver:  CFL=.5,  initial  density=l.l  Monolithic  Solver:  CFL  =  .5,  initial  density  =  1.1 


(a)  (b) 


Figure  17:  Numerical  profiles  for  bubble  volume  over  time  for  the  oscillating  bubble  problem  in  two  spatial  dimensions  with 
equation  of  state  p  =  Bp  generated  by  (a)  the  method  of  [5],  and  (b)  the  monolithic  solver. 

Next  consider  the  oscillating  bubble  problems  introduced  in  Section  2.2.  Here  we  use  the  two-way 
coupled  simulation  techniques  proposed  in  this  section  which  couple  the  incompressible  flow  solver  to  a 
full  compressible  flow  solver  that  includes  shocks  and  rarefactions,  albeit  a  somewhat  simplified  isothermal 
p  =  Bp  equation  of  state.  Figure  16(a)  shows  the  numerical  profiles  for  the  bubble  volume  over  time 
generated  using  the  partitioned  method  of  [5],  while  Figure  16(b)  shows  the  profiles  generated  using  our 
proposed  monolithic  solver  in  one  spatial  dimension.  Note  that  the  results  converge  to  the  “exact”  solution 
under  grid  refinement  in  both  cases.  Figure  17  shows  the  results  for  a  two  dimensional  oscillating  bubble. 
Note  that,  unlike  the  one  dimensional  case,  the  method  of  [5]  also  damps  the  solution  on  coarser  grids, 
although  converging  to  the  “exact”  solution  under  grid  refinement. 


4.  Constant  density  and  constant  pressure 

In  this  section  continuing  on  from  the  isothermal  assumption  from  Section  3.6,  we  additionally  make 
constant  density  and  constant  pressure  assumptions  to  arrive  at  our  final  method  for  simulating  bubbles. 

4-1.  Constant  density 

We  achieve  constant  density  by  redistributing  the  density  in  each  bubble  as  the  average  density  per 
bubble  before  the  implicit  pressure  solve,  exactly  as  is  done  for  the  partitioned  solver  in  Section  2.3.  Note 
that  a  spatially  constant  density  field  does  not  imply  that  pressure  inside  the  bubbles  is  spatially  constant  as 
well.  Figure  22  shows  the  numerical  profiles  for  the  bubble  volume  over  time  with  increasing  grid  resolutions 
for  the  oscillating  bubble  problems.  Note  that  the  profiles  converge  to  the  “exact”  solutions  under  grid 
refinement. 

4-2.  Constant  pressure 

We  now  present  a  simplified  monolithic  solver  which  solves  for  a  constant  pressure  pn+1  inside  the  bubble 
with  a  single  degree  of  freedom,  i.e.,  spatially  constant  but  time- varying.  In  matrix  terms,  this  corresponds 
to  taking  the  Poisson  matrix  for  the  full  system  and  collapsing  all  the  rows  and  columns  corresponding  to 
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Moving  Heavy  Incompressible  Droplet:  Density 


Moving  Heavy  Incompressible  Droplet:  Velocity 


(a)  density 


(b)  velocity 


X 

(c)  pressure 


Figure  18:  Numerical  results  for  the  moving  incompressible  droplet  example,  where  a  droplet  of  density  1000  kg/m3  is  travelling 
to  the  right  in  an  initially  stationary  compressible  fluid  with  equation  of  state  p  =  Bp  at  t  =  7.5  X  10— 4  s. 
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Moving  Light  Incompressible  Droplet:  Density  Moving  Light  Incompressible  Droplet:  Velocity 


x  x 

(a)  density  (b)  velocity 


Moving  Light  Incompressible  Droplet:  Pressure 


x 

(c)  pressure 


Figure  19:  Numerical  results  for  the  moving  incompressible  droplet  example,  where  a  droplet  of  density  10  kg/m3  is  travelling 
to  the  right  in  an  initially  stationary  compressible  fluid  with  equation  of  state  p  =  Bp  at  t  =  7.5  X  10— 4  s. 
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Shock  Impinging  on  Heavy  Incompressible  Droplet:  Density  Shock  Impinging  on  Heavy  Incompressible  Droplet:  Velocity 


(a)  density  (b)  velocity 


Shock  Impinging  on  Heavy  Incompressible  Droplet:  Pressure 


x 

(c)  pressure 


Figure  20:  Numerical  results  for  the  shock  impinging  on  a  heavy  incompressible  droplet  example  with  equation  of  state  p  =  Bp 
at  t  =  1.75  X  10— 3  s,  where  a  shock  wave  initially  located  at  x  =  .1  m  travels  to  the  right  impinging  on  an  incompressible 
droplet  of  density  1000  kg/m3  generating  both  reflected  and  transmitted  waves. 
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Shock  Impinging  on  Light  Incompressible  Droplet:  Density  Shock  Impinging  on  Light  Incompressible  Droplet:  Velocity 


(a)  density 


(b)  velocity 


Shock  Impinging  on  Light  Incompressible  Droplet:  Pressure 


x 

(c)  pressure 


Figure  21:  Numerical  results  for  the  shock  impinging  on  a  heavy  incompressible  droplet  example  with  equation  of  state  p  =  Bp 
at  t  =  1.75  X  lO-3  s,  where  a  shock  wave  initially  located  at  x  =  .1  m  travels  to  the  right  impinging  on  an  incompressible 
droplet  of  density  10  kg/m3  generating  both  reflected  and  transmitted  waves. 
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Constant  Density  Solver:  CFL  =  .125,  initial  density  =  1.1 


Constant  Density  Solver:  CFL  =  .125,  initial  density  =  1.1 


time 


(a) 


(b) 


Figure  22:  Numerical  profiles  for  bubble  volume  over  time  under  grid  refinement  generated  using  the  monolithic  solver  with 
density  redistribution  to  obtain  constant  density  as  well  as  an  equation  of  state  p  =  Bp  for  the  oscillating  bubble  problem  in 
(a)  one  spatial  dimension,  and  (b)  two  spatial  dimensions. 


cells  in  the  same  bubble  into  a  single  row  and  column,  adding  overlapping  matrix  elements.  This  can  also 
be  seen  as  summing  equation  (55)  over  the  entire  bubble  to  obtain 


w  — 1 - »n+1 

^  A  t2Bpn+^ 


EV'V 

ien 


(56) 


The  first  term  on  each  side  of  the  equality  simply  sums  over  the  number  of  cells  N  inside  the  bubble.  The 
last  term  can  be  modified  by  multiplying  and  dividing  by  the  volume  of  a  cell  Vc  and  converting  the  volume 
sum  to  a  surface  sum  along  the  MAC  grid  cell  faces  that  border  the  bubble.  If  the  average  normal  velocity 
on  all  these  faces  is  u  and  the  perimeter  of  all  these  faces  is  V,  we  obtain 


N 

At2Bpn+1 


nn+l 


TV  _  fbP 
At  ~  ~Vc 


(57) 


Most  of  the  terms  in  the  final  summation  vanish  since  Vp  is  zero  within  the  bubble,  leaving  only  terms 
corresponding  to  MAC  grid  faces  that  surround  the  bubble  volume.  For  each  of  these  faces  pn+1  equals  p 
as  defined  in  equation  (41).  Note  that  the  first  term  on  each  side  of  equation  (57)  is  based  on  the  number 
of  cells  within  the  bubble  and  the  second  term  on  each  side  of  the  equality  is  based  on  the  number  of  MAC 
grid  faces  surrounding  the  bubble.  Reminder,  pn+1  in  the  first  term  denotes  the  bubble  density. 

Note  that  the  constant  pressure  formulation  does  not  produce  pressure  gradients  within  the  bubble. 
Thus,  the  air  velocities  are  updated  through  a  second  projection  solve  as  per  equations  (15)  and  (16)  using 
the  boundary  incompressible  velocities,  as  outlined  in  Section  2.3.  Figures  23(a)  and  (b)  show  the  numerical 
profiles  for  the  bubble  volume  over  time  for  the  oscillating  bubble  problems  in  one  and  two  spatial  dimensions 
respectively.  Note  that  the  profiles  closely  match  those  shown  in  Figures  6  and  22,  verifying  the  correctness 
of  the  solver.  We  also  consider  the  case  where  the  initial  density  inside  the  bubble  is  1  kg/m3,  and  the  outside 
air  pressure  is  time-varying  as  pa tm(t)  =  Bf(t ),  where  /(£)  =  1  —  .2sin(27r£  +  7r/2).  We  refer  to  this  problem 
as  the  external  pressure  field  problem.  Figures  23(c)  and  (d)  show  the  resulting  numerical  profiles  for  the 
bubble  volume  over  time  generated  using  our  method  in  one  and  two  spatial  dimensions  respectively,  which 
converge  to  the  “exact”  solutions  under  grid  refinement.  Note  that  the  “exact”  solutions  are  computed  using 
equations  (63),  (70)  and  (77)  in  the  appendix  which  are  also  valid  for  time- varying  external  pressures. 

For  the  oscillating  bubble  and  external  pressure  field  problems,  we  also  considered  the  effects  of  surface 
tension  and  viscosity  for  a  range  of  parameters  from  high  to  low,  both  individually  and  in  combination. 


31 


a; 

£ 

o 


0.245 

0.24 

0.235 

0.23 

0.225 

0.22 

0.215 

0.21 

0.205 

0.2 


Constant  Pressure  Solver:  CFL  =  .125,  initial  density  =  1.1 
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(a) 


Constant  Pressure  Solver:  CFL  =  .125,  initial  density  =  1 


time 


(c) 


Constant  Pressure  Solver:  CFL  =  .125,  initial  density  =  1.1 


time 


(b) 


Constant  Pressure  Solver:  CFL  =  .125,  initial  density  =  1 


time 

(d) 


Figure  23:  Numerical  profiles  for  the  bubble  volume  over  time  under  grid  refinement  generated  by  the  constant  pressure 
formulation  where  all  pressure  unknowns  are  collapsed  to  a  single  degree  of  freedom  for  each  distinct  bubble  for  the  (a) 
oscillating  bubble  problem  in  one  spatial  dimension,  (b)  oscillating  bubble  problem  in  two  spatial  dimensions,  (c)  external 
pressure  field  problem  in  one  spatial  dimension,  and  (d)  external  pressure  field  problem  in  two  spatial  dimensions. 


Figures  24(a)  and  (b)  show  the  numerical  profiles  for  the  bubble  volume  over  time  for  both  these  problems 
in  two  spatial  dimensions  with  coefficients  a  =  .0728  kg/s2  and  p  =  .2  kg/ms.  To  heighten  the  effects  of 
surface  tension,  the  computational  domain  has  been  uniformly  scaled  down  by  10-3  m. 

As  further  validation  of  our  method,  we  also  considered  the  rising  bubble  examples  from  [18].  Consider  a 
computational  domain  of  [— 1  m,  1  m]  x  [— 1  m,  2  m]  which  is  initially  filled  with  water  except  for  a  circular 
air  bubble  of  radius  1/3  m  centered  at  the  origin  with  density  p  =  1.226  kg/m3.  The  effects  of  surface 
tension  and  viscosity  are  present  with  coefficients  p  =  .001137  kg/ms  and  a  =  .0728  kg/s2.  The  edges  of  the 
computational  domain  have  solid  wall  boundary  conditions.  Figure  25  shows  the  positions  of  the  air  bubble 
at  t  =  0,  t  =  .2,  t  =  .35  and  t  =  .5  seconds  under  grid  refinement.  Note  that  the  results  are  similar  to  those 
shown  in  Figures  3  and  4  from  [18].  The  finer  grid  computations  at  t  =  .35  s  and  t  =  .5  s  show  signs  of 
Kelvin- Helmholtz  instability,  as  noted  in  [18].  To  demonstrate  the  effects  of  surface  tension,  we  reduced  the 
computational  domain  to  [—.01  m,  .01  m]  x  [—.01  m,  .02  m]  and  the  radius  of  the  air  bubble  to  1/300  m. 
Figure  26  shows  the  positions  of  the  air  bubble  at  £  =  0,  £  =  .02,  t  =  .035  and  t  =  .05  seconds  under  grid 
refinement.  Note  that  the  results  are  similar  to  those  shown  in  Figures  1  and  2  from  [18]. 
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Constant  Pressure  Solver:  CFL  =  .5,  initial  density  =  1.1  Constant  Pressure  Solver:  CFL  =  .5,  initial  density  =  1 


(a)  (b) 


Figure  24:  Numerical  profiles  for  the  bubble  volume  over  time  under  grid  refinement  generated  by  the  constant  pressure  solver 
under  the  effects  of  surface  tension  and  viscosity  for  the  (a)  oscillating  bubble  problem  in  two  spatial  dimensions,  and  (b) 
external  pressure  field  problem  in  two  spatial  dimensions. 


5.  Momentum  conservation 

Consider  an  isolated  incompressible  droplet  with  no  ambient  pressure  forces.  When  solving  equation 
(4)  to  make  the  velocities  divergence  free,  a  constant  pressure  along  the  boundary  leads  to  a  net  force  of 
zero  implying  momentum  conservation  for  the  droplet.  When  updating  the  velocity  degrees  of  freedom  on 
regular  MAC  grid  faces  via  equation  (5),  the  interior  cell-centered  pressures  are  applied  in  a  conservative 
fashion  to  velocities  at  faces  bordering  interior  cells,  and  as  there  is  no  net  pressure  force  along  the  boundary 
the  projection  step  conserves  momentum  in  each  Cartesian  direction  independently.  That  is,  momentum 
is  conserved  during  the  projection  step  for  the  velocity  degrees  of  freedom  that  surround  the  cell-centered 
pressure  degrees  of  freedom. 

During  the  velocity  extrapolation  step,  if  one  computes  0- values  at  faces  by  averaging  cell-centered  <p- 
values,  then  some  of  the  faces  involved  in  the  momentum  conserving  projection  step  above  may  be  deemed 
outside  the  droplet  and  overwritten.  This  violates  momentum  conservation  and  can  be  seen  as  moving 
the  boundary  inward  by  one  grid  cell  replacing  the  proper  exterior  constant  pressure  Dirichlet  boundary 
condition  with  a  spurious  internal  pressure.  Thus  we  do  not  use  face- averaged  <p- values,  instead  labeling 
every  face  adjacent  to  an  interior  cell  center  as  interior  to  the  droplet.  Similarly,  the  viscous  solver  uses  the 
same  velocity  degrees  of  freedom  with  exterior  Neumann  boundary  conditions.  Currently,  we  do  not  use 
momentum  conserving  advection  [21]  as  the  largest  errors  for  momentum  conservation  occur  when  applying 
surface  tension. 

In  the  presence  of  surface  tension,  there  is  additional  momentum  added  for  any  velocity  degree  of  freedom 
that  lies  between  an  interior  and  an  exterior  cell.  We  compute  curvature  at  a  face  as  a  ^-weighted  average 
of  its  cell-centered  values,  i.e., 


„  _  Ki\4>i+1\  +  K*+l|0i| 

—  I  /  I  I  I  1  I  (ho) 

m\  +  Iw+il 

The  surface  tension  force  per  unit  area  at  this  face  is  sign(±l )ots,  where  sign(±l)  is  chosen  consistent  with 
the  outward  unit  normal  in  the  Cartesian  grid  directions  based  on  which  cells  are  interior  and  exterior  to  the 
level  set.  To  properly  conserve  momentum,  the  net  force  due  to  surface  tension  should  sum  to  zero  for  each 
independent  bubble  and  droplet  independently  along  each  of  the  Cartesian  directions.  To  enforce  this  for 
each  independent  bubble  and  droplet,  we  compute  the  total  surface  tension  force  per  unit  area  along  each 
Cartesian  direction  as  <7fttotal  =  sign(d=l)crtt/,  and  subtract  it  off  from  the  corresponding  force  for  each 
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Figure  25:  A  large  air  bubble  rising  in  water  with  solid  wall  boundary  conditions  under  grid  refinement.  Note  that  the 
computations  at  time  t  =  .5  s  show  signs  of  Kelvin-Helmholtz  instability. 
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Figure  26:  Level  set  profiles  generated  using  the  constant  pressure  solver  for  a  small  air  bubble  rising  in  water  with  solid  wall 
boundary  conditions  under  the  effects  of  surface  tension  and  viscosity. 
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face  in  a  curvature- weighted  fashion.  That  is,  the  new  jump  for  each  face  becomes 


crftjew  =  sign(±l)crK/  —  a/^total  (59) 

2^f  \Kf\ 

Obviously  alternatives  to  equation  (59)  exist,  and  it  is  not  the  form  of  the  correction  but  the  fact  that  a 
correction  needs  to  be  made  in  order  to  conserve  momentum  that  we  stress.  In  our  simulations  we  noticed 
that  the  net  surface  tension  force  along  the  boundary  is  close  to  zero  for  large  well-resolved  bubbles  and 
water  droplets.  However,  for  under-resolved  droplets  which  are  a  few  grid  cells  wide,  we  noticed  that  the  net 
surface  tension  force  can  be  so  far  from  zero  that  droplets  can  even  change  directions  in  mid-air  violating 
conservation  of  momentum. 

Although  we  only  apply  the  surface  tension  correction  to  closed  surfaces,  it  appears  that  one  can  use  a 
similar  strategy  for  open  surfaces  connected  to  boundaries  as  well.  Consider  for  example  a  single  circular 
bubble,  if  this  bubble  is  cut  in  half  then  it  follows  from  the  conservation  of  momentum  that  the  force  on  the 
top  half  of  the  bubble  must  equal  the  force  on  the  bottom  half.  However,  if  we  redrew  the  bottom  half  of 
the  bubble  to  any  arbitrary  curve,  then  the  net  force  on  this  new  bottom  half  must  still  be  equal  to  the  net 
force  on  the  top  half.  Essentially,  the  net  force  on  the  bottom  half  of  the  bubble  is  independent  of  its  shape 
and  is  somehow  tied  to  the  boundary  conditions  on  the  top  half.  Therefore,  for  an  open  surface  connected 
to  the  boundary  of  the  domain,  with  additional  knowledge  about  contact  angles,  one  could  likely  compute  a 
consistent  measure  for  the  net  force  on  the  open  surface. 


6.  Complex  bubble  breakup 

To  summarize,  our  method  for  simulating  bubbles  in  free  surface  incompressible  flows  is  as  follows.  We 
use  the  level  set  method  [24]  to  track  the  interface  between  the  bubbles  and  the  incompressible  fluid.  Initially, 
each  bubble  is  assigned  a  density  (or  mass)  which  is  advected  using  the  unconditionally  stable  conservative 
semi-Lagrangian  scheme  of  [21].  We  advect  the  bubble  density  using  air  velocities  which  are  constructed 
and  maintained  in  a  separate  velocity  field,  although  one  could  also  use  velocities  extrapolated  from  the 
incompressible  flow  for  increased  efficiency  but  lower  accuracy.  The  velocity  used  for  level  set  advection  is  a 
hybrid  between  the  incompressible  flow  velocity  and  the  air  velocity  field.  Despite  more  accurate  velocities, 
numerical  smearing  and  other  errors  cause  the  location  of  the  zero  level  set  and  the  location  of  the  non-zero 
bubble  densities  to  drift  apart  over  time.  We  address  this  issue  as  follows.  First  we  compute  the  total 
mass  that  belongs  to  a  bubble  as  the  sum  of  all  the  mass  inside  the  bubble  and  all  the  mass  closest  to 
that  bubble.  Then  we  use  a  flood  fill  algorithm  on  that  bubble  to  identify  all  grid  cells  belonging  to  that 
connected  component.  The  volume  of  this  connected  component  is  carefully  computed  using  a  piecewise 
linear  reconstruction  of  the  level  set  as  outlined  in  [20] .  The  mass  is  then  uniformly  redistributed  inside  the 
bubble  to  obtain  a  spatially  constant  bubble  density.  Before  advecting  the  incompressible  velocities  and  the 
air  velocities,  they  are  extrapolated  across  the  interface  in  order  to  define  ghost  node  values  using  the  fast 
extension  method  of  [1].  The  two  velocity  fields  are  then  independently  advected  using  the  second  order 
accurate  MacCormack  method  [27].  The  level  set  function  is  advected  using  the  particle  level  set  method 
of  [7]  and  the  semi-Lagrangian  advection  scheme  of  [8].  To  keep  the  level  set  a  signed  distance  function  we 
use  the  modified  fast  marching  method  of  [23].  Note  that  we  also  compute  the  advection  terms  in  a  thin 
band  of  ghost  cells  so  they  are  adequately  defined  when  the  interface  moves.  Also  note  that  we  passively 
advect  the  bubble  mass  that  is  not  close  to  any  bubble  using  the  incompressible  flow  velocities  in  order  to 
accurately  track  sub-grid  level  details.  When  a  bubble  is  about  to  merge  with  the  ambient  air  near  the  free 
surface,  dynamic  events  can  cause  the  bubble  to  open  in  one  time  step  and  close  in  the  next.  To  robustly 
track  the  bubble  density  in  these  cases,  we  keep  advecting  the  bubble  mass  even  if  the  bubble  merges  with 
the  ambient  air  near  the  free  surface. 

Viscosity  in  the  water  is  treated  implicitly  using  Neumann  boundary  conditions  at  the  interface  that  the 
derivative  of  each  component  of  the  incompressible  velocity  is  zero.  As  described  in  Section  2.1,  this  is  based 
on  the  assumption  that  the  dynamics  inside  the  air  bubbles  contain  little  momentum  and  hence,  they  cannot 
absorb  any  viscous  momentum  from  the  liquid.  Thus,  the  jump  in  pressure  due  to  viscosity  is  also  zero  since 


36 


the  normal  component  of  the  viscous  stress  vanishes  across  the  interface.  Finally,  as  we  assume  constant 
viscosity  in  the  incompressible  fluid,  the  equations  for  the  different  components  of  the  incompressible  velocity 
decouple  as  in  equations  (6)  and  (7)  for  two  spatial  dimensions.  Note  that  the  level  set  is  advected  to  its  new 
time  £n+1  position  before  the  viscous  update  due  to  the  need  to  prescribe  interface  boundary  conditions. 


(a)  (b) 


Figure  27:  Velocity  field  at  t  =  .5  s  on  a  80  X  120  grid  for  the  rising  bubble  example  from  [18]  (and  Figure  25)  when  the  air 
velocities  are  (a)  wiped  out  at  the  end  of  every  time  step  and  computed  from  the  boundary  incompressible  velocities  using  a 
second  projection  solve,  and  (b)  advected  forward  in  time  and  updated  using  pressure  gradients  from  a  second  projection  solve. 
The  incompressible  velocities  are  shown  in  blue  and  the  air  velocities  are  shown  in  red.  Note  that  the  velocity  field  in  (b) 
appears  much  more  continuous  and  smooth  compared  to  the  velocity  field  in  (a). 


For  the  implicit  step  of  our  method,  we  use  a  coupled  solve  between  a  single  degree  of  freedom  pressure 
for  the  bubble  and  the  surrounding  incompressible  flow  by  solving  equation  (57),  as  described  in  Section  4.2. 
The  ambient  air  is  taken  to  be  a  Dirichlet  boundary  condition  of  patm  =  101,325  Pa  for  which  we  use  the 
second  order  cut-cell  method  of  [11].  In  the  presence  of  surface  tension,  the  appropriate  cr  hi  jumps  are  added 
to  pressure  values  near  the  interface.  These  pressure  jumps  are  carefully  computed  noting  that  the  net 
surface  tension  force  on  each  bubble  and  each  water  droplet  must  be  zero,  as  described  in  Section  5.  We  use 
this  pressure  to  update  the  incompressible  velocities  via  equation  (5).  This  provides  a  very  stable  monolithic 
coupling  for  interactions  between  bubbles  and  the  surrounding  incompressible  flow.  The  air  velocities  inside 
the  bubbles  are  computed  from  the  boundary  incompressible  flow  velocities  using  a  second  projection  as 
in  equations  (15)  and  (16).  One  might  consider  wiping  out  the  air  velocities  at  the  end  of  every  time  step 
assuming  that  the  bubble  has  little  momentum  and  thus,  a  very  small  influence  on  the  incompressible  velocity 
field.  However,  this  destroys  the  temporal  continuity  of  the  overall  flow  field  as  shown  in  Figures  27(a)  and 
(b)  for  the  rising  bubble  example  from  Figure  25.  Note  that  all  the  vorticity  in  the  air  velocities  is  confined 
near  the  boundaries  in  Figure  27(a)  which  would  give  rise  to  a  noisy  level  set  as  the  simulation  progresses. 
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In  contrast,  Figure  27(b)  has  a  much  better  behaved  velocity  field  which  maintains  a  smooth  level  set  over 
time. 

Finally,  the  size  of  the  time  step  is  computed  using  equation  (10)  where  we  add  an  estimate  for  |  Vp|  using 
equation  (14)  which  accounts  for  the  change  in  velocity  over  a  time  step.  This  term  prevents  the  time  step 
from  becoming  excessively  large  when  velocities  are  near  zero  and  is  similar  in  spirit  to  the  idea  proposed 
in  [18]  for  body  forces  and  [19]  for  compressible  flow. 

6.1.  Numerical  results 

Consider  a  computational  domain  of  [— 1  m,  1  m]  x  [— 1  m,  2  m]  which  is  initially  filled  with  water  with 
density  1000  kg/m3  where  the  free  surface  is  located  at  y  =  1.5  m.  A  circular  air  bubble  of  radius  1/3  m 
is  centered  at  the  origin  with  density  1.364  kg/m3.  The  effects  of  surface  tension  and  viscosity  are  absent. 
Figure  28  shows  the  level  set  at  t  =  0,  t  =  .45,  t  =  .9,  t  =  1.2,  t  =  1.5,  t  =  1.8,  t  =  2.4,  t  =  2.7  and  t  =  3 
seconds  for  two  grids  of  resolutions  320  x  480  and  640  x  960  respectively.  Note  the  small  scale  details  that 
our  solver  is  able  to  resolve  and  accurately  track  over  time.  In  order  to  show  convergence  of  our  solver  under 
grid  refinement,  we  reduced  the  computational  domain  to  [—.01  m,  .01  m]  x  [—.01  m,  .02  m]  and  the  radius 
of  the  bubble  to  1/300  m,  and  included  the  effects  of  surface  tension  and  viscosity  with  coefficients  a  =  .0728 
kg/s2  and  /a  =  .001137  kg/ms.  The  initial  density  of  the  air  bubble  is  1.227  kg/m3.  Figure  29  shows  the 
level  set  at  t  =  0,  t  =  .02,  t  =  .04  and  t  =  .08  seconds. 

As  further  validation,  we  also  simulated  the  rising  bubble  experiment  corresponding  to  Fig  1A  from  [16] 
where  a  bubble  of  radius  .0061  m  and  initial  density  1  kg/m3  rises  in  a  liquid  of  density  875.5  kg/m3.  The 
computational  domain  is  [—.1464  m,  .1464  m]  x  [—.0732  m,  .366  m]  and  the  bubble  is  initially  located  at 
the  origin.  The  coefficient  of  viscosity  is  /r  =  .118  kg/ms  and  the  coefficient  of  surface  tension  is  a  =  .0322 
kg/s2.  The  edges  of  the  computational  domain  have  slip  solid  wall  boundary  conditions.  This  example  was 
also  simulated  in  [31,  15]  and  our  solver  gives  similar  results.  Figure  30(a)  shows  the  bubble  at  time  t  =  .6 
s  on  a  512  x  768  grid  along  with  the  streamlines  both  inside  the  bubble  and  in  its  wake.  Figure  30(b)  shows 
the  time  evolution  of  the  position  of  the  center  of  mass  of  the  bubble.  We  compared  this  data  with  the  linear 
best  fit  for  .6  <  t  <  .8  seconds.  The  slope  of  the  linear  best  fit  is  .2114  m/s  whereas  the  expected  slope  is 
.215  m/s.  We  note  that  the  streamlines  do  not  match  the  experimentally  observed  values  due  to  the  fact 
that  a  highly  simplified  single  pressure  degree  of  freedom  model  is  used  for  the  air,  however,  the  steady  state 
bubble  rise  speed  is  close  to  the  experimentally  measured  value  in  [16]  and  the  steady  state  bubble  shape  is 
similar  to  that  observed  in  [16]  and  computed  in  [31,  15]. 

6.1.1.  Object  interaction 

When  updating  the  incompressible  flow  velocities  during  advection,  in  the  presence  of  objects,  we  set 
velocity  Dirichlet  boundary  conditions  at  cell  faces  whose  centers  he  inside  an  object  as  these  cell  faces  have 
well-defined  velocities  determined  by  the  object  velocity.  We  then  advect  every  face  obtaining  a  well-defined 
incompressible  flow  velocity  field  at  time  £*.  As  we  consider  interaction  with  objects  that  are  fully  submerged 
initially,  we  initialize  j)  inside  these  objects  to  be  the  value  when  no  objects  were  present,  i.e.,  treat  them  as 
water.  Subsequently  when  advecting  </>,  we  also  update  j)  at  grid  cells  whose  centers  he  inside  an  object.  Air 
velocities  are  handled  similar  to  the  incompressible  flow  velocities  during  advection  with  Dirichlet  boundary 
conditions  at  cell  faces  inside  objects.  As  described  in  [21]  (see  also  [14]),  for  advecting  the  air  density, 
we  modify  the  forward  and  backward  ray  casting  to  stop  when  it  hits  an  object.  Interpolation  weights 
are  computed  at  the  surface  point,  where  weights  coming  from  cells  inside  an  object  are  discarded  and  the 
remaining  weights  are  rescaled  to  sum  to  1.  During  the  viscous  solve,  we  set  Dirichlet  boundary  conditions  at 
cell  centers  that  he  inside  objects.  When  solving  for  the  pressure  to  make  the  incompressible  flow  velocities 
divergence  free,  we  set  Neumann  boundary  conditions  at  cell  faces  whose  centers  he  inside  objects.  Similarly, 
during  the  second  projection  step  for  updating  the  air  velocities,  we  set  Neumann  boundary  conditions  at 
cell  faces  whose  centers  he  inside  objects. 

Consider  a  computational  domain  of  [— 1  m,  1  m]  x  [— 1  m,  2  m]  which  is  initially  filled  with  water  with 
density  1000  kg/m3  where  the  free  surface  is  located  at  y  =  1.5  m.  A  circular  air  bubble  of  radius  1/3  m 
is  centered  at  the  origin  with  density  1.364  kg/m3.  The  domain  has  nine  circular  objects  with  four  objects 


38 


t  =  2.4  s 


t  =  2.7  s 


tm  3  s 


Figure  28:  An  air  bubble  rising  in  a  water  column  with  a  free  surface.  The  effects  of  surface  tension  and  viscosity  are  absent. 
Note  the  small  scale  details  that  our  solver  is  able  to  resolve  and  accurately  track  over  time. 
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Figure  29:  Level  set  profiles  under  grid  refinement  for  an  air  bubble  rising  in  a  water 
computational  domain.  The  effects  of  surface  tension  and  viscosity  are  present. 


column  with  a  free  surface  in  a  small 
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Figure  30:  Rising  bubble  example  from  Fig.  1A  of  [16],  (a)  computed  steady  state  bubble  shape  and  the  streamlines  inside  the 
bubble  and  those  in  its  wake,  and  (b)  time  evolution  of  the  position  of  the  center  of  mass  of  the  bubble. 


located  at  (—.6  +  .4fc,  .5),  where  k  =  {0, 1,  2,  3}  and  the  other  five  objects  located  at  (—.8  +  .4fc,  .9),  where 
k  =  {0,1,  2,  3, 4}.  Figure  31  shows  the  level  set  at  t  =  0,  t  =  .25,  t  =  .5,  t  =  .75,  t  =  1,  t  =  1.25,  t  =  1.5, 
t  =  1.75  and  t  —  2.1  seconds  for  two  grids  of  resolutions  320  x  480  and  640  x  960  respectively.  As  can  be 
seen,  the  objects  break  up  the  bubble  into  a  large  number  of  small  bubbles  which  our  solver  is  able  to  resolve 
and  efficiently  track  over  time.  We  also  show  an  example  where  the  larger  bubbles  deform  while  the  smaller 
bubbles  are  able  to  preserve  their  shape  due  to  surface  tension  forces  by  reducing  the  computational  domain 
to  [—.1  m,  .1  m]  x  [—.1  m,  .2  m],  the  radius  of  the  bubble  to  1/30  m  with  initial  density  1.24  kg/m3,  and 
adding  the  effects  of  surface  tension  and  viscosity  with  coefficients  a  =  .0728  kg/s2  and  /i  =  .001137  kg/ms. 
Figure  32  shows  the  level  set  at  t  =  0,  t  m  #1,  t  =  .2,  t  =  .3,  t  =  .4,  t  =  .5,  t  =  .6,  t  =  .7  and  t  =  .8  seconds 
for  two  grids  of  resolutions  160  x  240  and  320  x  480  respectively.  Note  that  the  smaller  bubbles  remain 
spherical  because  of  larger  surface  tension  forces  while  the  larger  bubbles  readily  deform.  In  order  to  show 
convergence  under  grid  refinement  in  the  presence  of  objects,  we  reduced  the  computational  domain  to  [—.02 
m,  .02  m]  x  [—.02  m,  .04  m]  and  the  radius  of  the  bubble  to  1/150  m,  and  reduced  the  number  of  objects 
to  five.  Two  objects  are  located  at  (—.007  +  .014 &,  .01),  where  k  =  {0, 1}  and  the  other  three  are  located  at 
(— .014  +  .014 &,  .018),  where  k  =  {0, 1,  2}.  The  effects  of  surface  tension  and  viscosity  are  present  with  the 
same  coefficient  values  as  above.  The  initial  density  of  the  air  bubble  is  1.229  kg/m3.  Figure  33  shows  the 
level  set  at  t  =  0,  t  =  .075,  t  =  .16  and  t  =  .225  seconds. 

6.1.2.  Three  spatial  dimensions 

Consider  a  computational  domain  of  [— 1  m,  1  m]  x  [— 1  m,  2  m]  x  [— 1  m,  1  m]  which  is  initially  filled 
with  water  with  density  1000  kg/m3  where  the  free  surface  is  located  at  y  =  1.5  m.  An  air  bubble  of  radius 
1/3  m  and  density  1.364  kg/m3  is  centered  at  the  origin.  The  effects  of  surface  tension  and  viscosity  are 
absent.  To  break  up  the  bubble  into  a  large  number  of  smaller  bubbles,  the  domain  also  has  25  spherical 
objects  of  radius  .15  m  centered  at  (—.8  +  .4i,  .5,  —.8  +  .4j),  where  i,  j  =  {0, 1,  2,  3, 4}.  Figure  34  shows  the 
level  set  at  t  =  0,  t  =  .7,  t  =  1  and  t  =  1.4  seconds  simulated  using  a  grid  of  resolution  256  x  384  x  256.  A 
close  up  of  the  level  set  at  t  =  1  seconds  is  shown  in  Figure  35  illustrating  the  large  amount  of  topological 
detail  that  our  solver  is  able  to  resolve  and  accurately  track  over  time. 
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Figure  31:  An  inviscid  air  bubble  rising  in  a  water  column  with  a  free  surface  in  the  presence  of  objects.  The  objects  break  up 
the  bubble  into  a  large  number  of  small  bubbles  which  our  solver  is  able  to  resolve  and  efficiently  track  over  time. 
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t  —  Os 


t  —  .  3  s  t  =  A  s 


Figure  32:  An  air  bubble  rising  in  a  water  column  with  a  free  surface  in  the  presence  of  objects  in  a  medium-sized  computational 
domain.  The  effects  of  surface  tension  and  viscosity  are  present.  Note  that  the  smaller  bubbles  remain  spherical  because  of 
larger  surface  tension  forces  while  the  larger  bubbles  readily  deform. 
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Figure  33:  Level  set  profiles  under  grid  refinement  for  an  air  bubble  rising  in  a  water  column  with  a  free  surface  in  the  presence 
of  objects  in  a  small  computational  domain.  The  effects  of  surface  tension  and  viscosity  are  present. 
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Figure  34:  An  air  bubble  rising  in  a  water  column  with  a  free  surface  in  the  presence  of  objects  in  three  spatial  dimensions. 
The  effects  of  surface  tension  and  viscosity  are  absent. 
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Figure  35:  A  close  up  of  the  three  dimensional  rising  bubbles  at  t  =  1  seconds  (see  Figure  34(c)).  Note  the  large  amount  of 
topological  detail  that  our  solver  is  able  to  resolve  and  accurately  track  over  time. 


7.  Conclusion 

We  designed  a  method  for  simulating  air  bubbles  in  free  surface  incompressible  flows.  To  formulate 
our  method,  we  first  proposed  a  straightforward  partitioned  solver  based  on  mass  tracking.  We  showed  that 
such  an  approach  suffers  from  stability  issues  which  have  characteristics  similar  to  partitioned  (as  opposed  to 
monolithic)  methods  for  solid-fluid  coupling  [13,  26,  12].  These  issues  can  be  alleviated  using  outer  iterations 
on  the  partitioned  solver,  although  the  computation  time  increases  drastically  because  each  time  step  can 
require  ten  or  more  Poisson  solves.  Hence,  we  took  a  monolithic  approach  for  the  air-water  problem  similar 
to  the  solid-fluid  coupling  in  [13,  26,  12]  as  motivated  by  [19].  To  design  this  approach,  we  revisited  the 
partitioned  solver  of  [5]  for  coupling  compressible  and  incompressible  flow  and  devised  a  monolithic  solver 
using  the  ideas  from  [19]  to  couple  together  incompressible  flow  with  fully  non-linear  compressible  flow 
including  shocks  and  rarefactions.  We  then  simplified  this  approach  greatly  to  make  this  approach  in  line 
with  our  straightforward  partitioned  approach  for  simulating  bubbles.  This  was  achieved  by  setting  both  the 
bubble  density  and  the  bubble  pressure  to  be  spatially  constant,  although  time- varying.  We  demonstrated 
the  accuracy  and  robustness  of  this  method  on  test  problems  as  well  as  more  realistic  problems  in  both  one, 
two  and  three  spatial  dimensions.  In  the  future,  we  would  like  to  couple  this  method  with  those  for  modeling 
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cavitation  [4]  to  simulate  real-world  problems  such  as  ship  propellers  or  breaking  waves. 
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Appendix 


Consider  an  infinitesimal  element  with  volume  dfl  and  force  per  unit  volume  Vp,  implying  that  the  total 
force  on  this  element  is  VpdU.  The  work  done  in  displacing  this  element  in  an  infinitesimal  time  interval  dt 
is  given  by  VpdU  •  dx  =  VpdU  •  vdt.  The  total  work  done  dW  in  the  time  interval  dt  is  the  integral  of  this 
work  over  the  entire  domain,  i.e. , 

dW  =  f  Vp  •  vdQdt  =  f  V  •  (pv)dQdt  (60) 

Jn  Jn 

since  V  •  v  =  0.  Using  the  divergence  theorem,  this  integral  is  equivalent  to  the  following  surface  integral 

dW  =  <J)  pv  •  ndSdt  (61) 

Jdn 

which  we  use  in  the  following  subsections. 

Appendix-I 

Consider  the  one  dimensional  oscillating  bubble  problem  introduced  in  Section  2.2  and  Figure  1.  Since 
the  system  is  symmetric  about  the  midpoint  we  only  consider  the  right-half  of  the  domain.  In  one  spatial 
dimension,  water  being  incompressible  has  a  spatially  constant  velocity.  Since  there  is  no  mass  transfer 
between  the  air  bubble  and  the  water  this  is  also  the  velocity  of  the  bubble-water  interface,  or  the  rate  of 
change  of  the  radius  of  the  bubble.  Thus  at  time  £,  v(t)  =  R{t)  at  the  interface  of  the  bubble  of  radius  R{t). 
If  l  is  the  length  of  the  water  region,  then  pjl  is  its  mass.  Let  Pb(t)  be  the  pressure  inside  the  bubble  and 
Patm (t)  be  the  pressure  in  the  air  at  time  t.  The  total  force  on  the  water  is  given  by  Pb(t)  —  pa tm(^)-  This 
must  equal  the  mass  times  the  acceleration  R{t)  of  water,  i.e., 

Pb(t)  -  patm(t)  =  pilR(t)  (62) 

Using  an  equation  of  state  p  =  Bp ,  it  follows  that  Pb(t)  =  Bpb(t )  =  BM/R{t ),  where  M  is  the  constant 
mass  of  the  bubble.  Substituting  this  in  equation  (62)  gives 

=  ~Patm{t))  (63) 

Appendix-11 

Consider  the  oscillating  bubble  problem  in  two  spatial  dimensions  as  shown  in  Figure  2.  Since  the  total 
volume  of  water  is  conserved,  the  radius  of  the  water  sphere  Rw(t )  is  dependent  on  the  radius  of  the  bubble 
R(t).  Let  R°  be  the  initial  radius  of  the  bubble  and  R ^  be  that  of  the  water  drop.  Then  conservation  of 
volume  yields 

n((R°w )2  -  OR0)2)  =  tt (Rw(t)2~R(t)2) 

Rw{t)  =  y/(iQ2  +  RW  ~  (-R0)2 

Rw(t )  =  \/a2  +  R(t )2,  where  a 2  =  ( R ^)2  —  ( R °)2  (64) 

In  fact,  the  volume  is  conserved  for  any  annulus  with  inner  radius  R(t)  and  outer  radius  r,  so  that  7r (r2—R(t)2) 
is  constant  and 

rr  =  R(t)R(t)  (65) 

At  any  point  in  time  the  kinetic  energy  of  an  annulus  of  infinitesimal  thickness  dr  and  radius  r  is  given  by 
|(p/27rrdr)  •  v(r)2,  where  v(r )  =  r  =  R{t)R{t)/r  from  equation  (65).  The  total  kinetic  energy  of  water  is 
obtained  by  integrating  this  from  the  radius  of  the  bubble  R{t)  to  the  outer  radius  of  water  Rw(t ),  i.e., 

Rw(t) 

K.E.{t)  =  pm  j  (. R(t)R{t)/rfrdr  =  itPl(R(t)R{t))2  In (Rw(t)/R(t))  (66) 

m 
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Equating  the  work  done  with  the  kinetic  energy,  and  noting  that  equation  (61)  applied  to  our  annulus  gives 
two  terms  of  the  form  2tt rr  results  in 

t 

j pb(T)2nR(T)R(T)d,T  -  patm('r)27ri?TO(r )Rw(t)<1t  =  7 xpi(R(t)R(t))2  In (Rw(t)/R(t))  (67) 

0 

From  equation  (65),  Rw(t)Rw(t)  =  R(t)R(t)  for  any  time  t  and  thus 


t 

2  J  (Pb(r)  —  Patm  (t))R(t)R(t)cIt  =  pi(R(t)R(t)f,\n.(y^/o?  +  Riff/R(t)) 


(68) 


Differentiating  both  sides  with  respect  to  t  and  simplifying  gives 

Pb(t)  Patm (t)  =  pi(R(t)R(t)  +  R{t)2)  ln(\/ a2  +  R(ff/R(t))  -  .5a2 piR(t)2 /(a2  +  R{t)2)  (69) 

Substituting  Pb(t)  =  Bpb(t)  =  BM/(itR(t)2)  and  rearranging  gives 

.5a2  \ 


piR(t)R(t )  In 


^a2  +  R(t)2 

m 


+  piR(t )  In 


x /a2  +  R(t )2 

m 


BM 


a 2  +  R(t)2  j  7r R(t)2 


=  0  (70) 


Appendix-Ill 

Finally,  consider  the  three  dimensional  case.  Equation  (64)  becomes 

Rw(t)  =  y/a3  +  it*(£)3,  where  a3  =  (i?3  )3  -  ( R °)3 


equation  (65)  becomes 


r2r  =  R(t)2R{t) 


and  equation  (66)  becomes 


(71) 

(72) 


Rw(t) 

K.E.(t)  =  2ttPi  J  (R(t)2R{t)/r2)2r2dr  =  2npIR(t)AR(t)2{l/R(t)-l/Rw(t))  (73) 

m 

Equating  the  work  done  with  the  kinetic  energy,  and  noting  that  equation  (61)  applied  to  our  thickened  shell 
gives  two  terms  of  the  form  4tt r2r  results  in 


i 

J p\,(t )4ttR(t )2 R(t )dr  -  p3tm{T)A-KRw(r)2 Rw(r)dT  =  2npIR(t)4R(tf  (1/ R{t)  -  l/Rw(t))  (74) 


From  equation  (72),  i?^(t)i7w(t)  =  R2(t)R(t)  for  any  time  t  and  thus 


t 

■  J (Pb(r)  -  patm(r))R(T)2 R(r)dr  =  PlR(t)4R(tf  (l /R{t)  -  1/ ^a3  +  i?(i)3)  (75) 


Differentiating  both  sides  with  respect  to  t  and  simplifying  gives 


Pb(t)  -Patm(t)  =  piR(t)  R(t)  - 


m 2 


ya3  +  Wf, 


+  PiR{t)2  o  + 


m2 


2  R(t) 


Substituting  Pb(t)  =  Bpb(t)  =  B M / (|7r R(t)3)  and  rearranging  gives 


piR{t)  R(t)  - 


R{t)2 


2/ya3+i?(i)3 


+  piR{t)2  -  + 


R{tf 


2  1  2/y«3  +  y^+w 

2  m  \ 


BM 


(76) 


2/y«3  +  W)5  y«3  +  rW  I  f  7T R(t)3 


+  Pam(t)  =  0 
(77) 
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